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Dynamic  simulation  has  become  a necessary  tool  in  the  design  and  analysis 
process  of  mechanical  systems.  This  study  is  concerned  with  improving  on  the  current 
technologies  associated  with  the  simulation  of  mechanical  systems.  One  issue  considered 
deals  with  the  modeling  of  flexible  components.  Dynamic  simulation  packages  often 
utilize  the  assumed  modes  approach  for  the  modeling  of  flexible  components.  A 
relatively  new  idea  in  reduced  component  modeling  is  to  generate  reduced  components 
based  on  the  projection  of  system  level  modes  onto  the  components.  This  approach 
yields  component  models  that  incorporate  effects  of  the  inboard  and  outboard  boundary 
conditions.  This  methodology  is  extended  to  use  Krylov  system  level  modes  in  lieu  of 
the  conventional  normal  modes  as  a projection  basis.  Numerical  results  indicate  that  the 
reduced  Krylov  system  provides  an  excellent  alternative  since  it  does  not  require  the 
solution  of  the  full-order  eigenvalue  problem. 
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This  work  also  introduces  a switched  system  to  achieve  order  reduction  for  linear 
structural  systems.  The  switched  system  model  is  derived  based  on  modal  participation 
factors  calculated  from  an  input  forcing  function  and  initial  conditions  at  each  switch 
step.  Stability  issues  for  switched  systems  are  considered  followed  by  the  necessary 
conditions  of  stability  for  the  proposed  switched  model.  Numerical  simulations 
demonstrate  that  an  error  indicator  based  on  modal  participation  factors  provides  an 
effective  switching  strategy. 

A second  issue  addressed  is  concerned  with  the  numerical  solution  procedures  for 
the  governing  equations  of  some  mechanical  systems.  Resulting  governing  equations  can 
be  defined  in  a differential  algebraic  form  and  their  solutions  are  often  subject  to 
constraint  violation  drift.  In  this  work,  constraint  violation  formulations  are  derived 
within  the  framework  of  nonlinear  control  theory.  Using  Lie  algebraic  methods,  it  is 
possible  to  interpret  many  constraint  stabilization  methods  as  nonlinear  controllers. 
Alternatively,  somewhat  simpler  constraint  stabilization  methods  can  be  derived  using 
energy  integrals  associated  with  the  dynamical  system.  These  interpretations  present  a 
framework  that  allows  for  the  development  of  alternative  stabilization  methods  that  are 
simple  and  effective.  A non-trivial  numerical  example  is  provided  that  considers 
constraint  stabilization  between  two  subassemblies:  a highly  maneuverable  multi- 
wheeled vehicle  (HMMWV)  and  a weapon  assembly. 
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CHAPTER  1 
INTRODUCTION 


1.1  Motivation 

Dynamic  simulation  of  mechanical  systems  is  becoming  increasingly  important  in 
the  area  of  engineering  design  and  analysis.  The  motivation  for  dynamic  simulation 
stems  from  many  aspects  of  the  engineering  industry.  As  industry  needs  become  more 
and  more  advanced,  so,  too,  do  the  technologies  required  to  meet  those  needs.  The 
design  and  cost  of  simple  mechanisms  can  be  greatly  improved  with  the  aid  of  dynamic 
or  virtual  simulation.  Tests  and  analyses  can  be  performed  to  aid  in  feasibility  studies 
and  cost  estimations  before  the  design  is  sent  to  production.  Dynamic  simulation  is  also  a 
useful  tool  in  the  design  of  more  complex  structures,  such  as  automobiles  and  spacecraft. 
Predictions  of  these  system  responses  to  forces  or  disturbances  can  be  used  to  streamline 
design  concepts  and  aid  in  determining  system  specifications. 

Mechanical  systems,  in  general,  can  be  modeled  as  multibody  systems.  That  is, 
the  systems’  models  are  made  of  a multiple  bodies  whose  relative  motions  are  described 
by  hinges  or  joints.  Depending  on  the  physical  system,  the  bodies  can  be  modeled  with 
or  without  flexible  effects.  More  specifically,  bodies  in  a multibody  system  are  modeled 
as  one  of  three  possible  body  types:  1)  rigid  bodies,  2)  flexible  bodies  that  undergo  small 
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elastic  (linear)  motions  and  3)  flexible  bodies  that  undergo  large  elastic  (nonlinear) 
motions.  In  this  dissertation,  only  the  first  and  second  body  types  are  considered. 

The  models  of  the  multibody  systems  described  above  contain  a coupled  set  of 
ordinary  differential  equations  (ODEs)  and  partial  differential  equations  (PDEs).  The 
ODEs  represent  the  rigid  body  motion  of  the  system  while  the  PDEs  represent  the 
flexible  motion.  For  the  flexible  bodies  that  undergo  small  motions,  the  PDEs  can  be 
approximated  using  a spatial  discretization  yielding  system  equations  described  solely  by 
ODEs.  The  connectivity  or  constraints  on  the  system’s  bodies  provides  additional 
equations  of  an  algebraic  form.  The  coupled  set  of  ODEs  and  algebraic  equations  results 
in  system  equations  described  by  a class  of  equations  called  differential-algebraic 
equations  (DAEs). 

Systems  requiring  a large  number  of  degrees-of- freedom  (DOFs)  result  in  DAEs 
of  extremely  large  dimensions.  The  cost  of  obtaining  numerical  solutions  of  these 
equations  is  related  to  a number  of  system  parameters  including  the  dimensionality  of  the 
system.  Order  reduction  of  these  systems  without  loss  of  simulation  accuracy  is  a non- 
trivial task  and  is  a major  focus  of  this  dissertation. 

A second  topic  of  this  dissertation  deals  with  the  numerical  solution  of  DAEs 
associated  with  multibody  dynamical  (MBD)  systems.  Unlike  numerical  solutions  to 
ODEs,  solutions  to  these  DAEs  are  quite  complex.  Due  to  the  nature  of  DAEs,  problems 
associated  with  constraint  violation  drift  can  arise.  The  error  or  drift  is  conceived  in  the 
formulation  of  the  DAE.  In  order  for  a MBD  system’s  equations  to  be  constructed  as 
DAEs,  the  constraint  accuracy  is  compromised  by  considering  the  constraint 


accelerations  instead  of  the  constraints  themselves.  Stabilization  of  these  constraints  is 
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the  central  topic  in  numerous  research  articles.  Recently,  it  was  shown  that  the 
stabilization  of  the  constraints  can  be  thought  of  as  a controls  problem  that  utilizes  partial 
feedback  linearization.  This  interpretation  of  constraint  stabilization  and  its  applications 
plays  an  important  role  in  this  dissertation. 

1.2  Objectives 

The  present  study  seeks  to  improve  on  the  standard  practices  for  the  dynamic 
simulation  of  mechanical  systems,  more  specifically,  mutliflexible  body  structural 
systems.  The  main  objectives  are  to  (1)  improve  on  existing  model  reduction  techniques 
and  develop  novel  reduction  techniques  based  on  model  switching  without  the  sacrifice  of 
simulation  accuracy  and  (2)  demonstrate  the  value  of  constraint  stabilization  via 
nonlinear  control  methodologies. 

In  Chapter  2,  a review  of  the  literature  used  in  the  development  of  these  objectives 
is  provided.  Topics  covered  include  general  model  reduction,  constraint  stabilization,  and 
switched  systems. 

In  Chapter  3,  the  fundamentals  of  system  level  model  reduction  are  presented.  A 
discussion  on  the  Projection  and  Assembly  methods  for  articulating  systems  is  given. 
Details  of  the  basis  choices  are  given,  followed  by  numerical  examples  that  verify  the 
various  techniques’  benefits. 

In  Chapter  4,  the  concepts  of  constraint  stabilization  are  presented.  First,  the 
numerical  problems  associated  with  the  corresponding  DAEs  are  derived.  Stabilization 
techniques  are  then  discussed  and  an  interpretation  of  the  stabilization  techniques  shows 
that  the  problem  can  be  construed  as  a partial  feedback  linear  controller.  The  stabilization 
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is  then  verified  in  a non-trivial  numerical  example  involving  a military  vehicle 
(HMMWV)  and  a weapon  system. 

Chapter  5 approaches  the  model  reduction  problem  with  a novel  twist.  As  a first 
step,  only  linear  structural  systems  are  considered.  The  idea  is  to  model  linear  systems  as 
switched  systems.  Definitions  and  preliminaries  of  switched  systems  are  initially 
provided  as  well  as  additional  necessary  system  requirements.  A proposed  switched 
system  model  is  then  derived  for  linear  structural  systems.  Next,  a definition  of  stability 
is  given  along  with  simple  examples  of  stable  and  unstable  switched  systems.  Finally,  a 
proof  of  stability  is  derived  for  systems  satisfying  certain  parameter  specifications. 

In  Chapter  6,  numerical  examples  provide  testbeds  for  exploring  fundamental 
concepts  and  demonstrating  the  potential  for  switched  systems  in  linear  structural 
analysis.  A simple  two  DOF  spring  damper  system  is  provided  to  demonstrate  some  of 
the  elementary  effects  of  switched  system  modeling  using  the  proposed  switched  model. 
This  is  followed  with  a more  complex  example  of  a 30  DOF  cantilevered  beam.  This 
example  is  intended  to  present  some  of  the  benefits  of  using  modal  switching  on  linear 
structural  systems. 

In  Chapter  7,  summaries  and  conclusions  of  the  topics  discussed  in  this  study  are 
presented  along  with  some  suggestions  for  future  work. 


CHAPTER  2 
LITERATURE  REVIEW 


2.1  Introduction 

The  objective  of  this  study  is  to  improve  the  modeling  techniques  of  mechanical 
or  multiflexible  body  systems.  Formulating  accurate  governing  equations  of  motion  and 
obtaining  robust  solutions  to  these  equations  is  the  ultimate  goal.  The  search  to  meet  this 
goal  requires  a vast  number  of  subtopics  in  the  area  of  dynamic  simulation.  This  chapter 
reviews  literature  associated  with  some  of  these  subtopics. 

In  many  mechanical  systems,  flexible  components  can  be  modeled  using  linear 
elastic  theory.  Substructuring  and  modal  selection  techniques  are  examples  of  concepts 
developed  to  enhance  flexible  component  modeling.  Section  2.2  reviews  these  and  other 
developments  in  the  area  of  flexible  component  modeling.  Additional  topics  in 
component  modeling  are  associated  with  the  use  of  system  level  modes.  Various 
projection  and  assembly  methods  have  been  developed  for  the  purpose  of  generating 
component  modes  from  the  system  level.  Their  usefulness  in  linear  and  articulating 
systems  is  reviewed  in  Section  2.3. 

Resulting  governing  equations  of  mulitflexible  body  systems  are  of  the 
differential  algebraic  form.  Obtaining  numerical  solutions  to  the  DAEs  poses  theoretical 
difficulties  in  the  form  of  constraint  violation  drift.  Section  2.4  reviews  constraint 
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implementation  methods  and  the  various  techniques  which  seek  to  minimize  constraint 
violations. 

Structural  dynamical  systems  often  comprise  flexible  components  whose 
configurations  can  depend  on  time.  Conventional  modal  selection  methods  are  not 
capable  of  being  time-dependent,  which  makes  optimizing  these  modal  selection 
procedures  difficult.  Similar  difficulties  arise  in  choosing  minimal  modal  bases  for  linear 
structural  systems  whose  applied  forcing  functions  vary  in  direction.  These  issues 
motivate  the  idea  for  mode  switching  or  modeling  dynamical  systems  as  switched 
systems.  Section  2.5  reviews  the  mathematical  developments  and  applications  of 
switched  systems  in  engineering. 

2.2  Linear  System  Reduction 

2.2. 1 Flexible  Multibody  Dynamic  Formulations 

Most  flexible  multibody  systems  can  be  thought  of  as  of  systems  whose 
components  undergo  large  rigid  body  motions  with  superimposed  small  elastic  motions. 
The  flexible  motion  is  often  described  from  a reference  frame  on  the  component  in 
motion.  This  technique,  termed  floating  reference  frame  [57],  has  many  advantages  and  is 
a widely  accepted  procedure  in  dynamic  simulation  codes.  One  advantage  is  that  the 
modeling  of  flexible  components  is  independent  of  the  dynamic  simulation.  This  means 
that  the  parameters  (modal  parameters)  used  to  describe  the  flexible  characteristics  can  be 
generated  off  line,  before  the  simulation. 
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Flexible  components  can  be  viewed  as  distributed  parameter  systems  or 
continuous  systems.  In  many  cases,  these  components  can  be  modeled  with  linear  PDEs. 
For  obvious  numerical  reasons,  the  inclusion  of  solutions  to  the  infinite  dimensional 
PDEs  with  standard  dynamic  simulation  codes  is  not  feasible.  Approximation  techniques 
are  then  utilized  to  model  the  PDEs.  Two  of  the  most  common  methods  for  the 
approximation,  assumed  modes  method  and  finite  element  method  (FEM)  [3 1 ] are  spatial 
discretizations  methods  which  convert  the  infinite  dimensional  PDEs  into  a finite  set  of 
ODEs  [34],  FEM  is  quite  useful  for  systems  with  complex  geometries,  but  it  often 
generates  system  models  of  extremely  large  dimensions.  The  assumed  modes  approach, 
however,  requires  an  intuitive  selection  of  basis  modes  such  that  the  DOFs’  motions  are 
well  described.  Dynamic  simulation  packages  such  as  DADS  [14],  TREETOPS  [29], 
DISCOS  [11]  and  ADAMS  [50]  allow  for  flexible  components  to  be  modeled  using  the 
assumed  modes  approach.  It  is  for  this  reason,  that  the  current  study  concentrates  on  the 
modeling  of  linear  flexible  components  using  the  assumed  modes  approach. 

The  assumed  modes  method  approximates  the  deflections  of  elastic  components 
using  a finite  series  of  space-dependent  functions  multiplied  by  time-dependent  functions 
[15,  34,  57].  The  time-dependent  functions  are  actually  amplification  factors  of  the 
space-dependent  functions  and  are  additional  DOFs  or  generalized  coordinates  in  the 
systems  governing  equations.  Some  formulations  replace  the  functions  with  vectors, 
yielding  a more  simplistic  approximation  of  the  flexible  component  [58].  In  DADS,  a 
FEM  of  the  flexible  component  is  passed  to  an  intermediate  processor  that  creates  the 
necessary  flexible  data  based  on  the  vector  formulation.  This  study  will,  therefore,  focus 
on  formulations  using  basis  vectors  instead  of  basis  functions. 
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2.2.2  Component  Mode  Synthesis 

For  components  that  are  modeled  with  thousands  of  DOFs,  it  is  essential  to 
choose  a minimal  basis  set  that  captures  the  important  characteristics  of  the  system. 
Numerous  articles  have  been  written  on  how  best  to  select  the  basis  vectors.  One  of  the 
most  common  techniques  for  modeling  linear  structures  is  component  mode  synthesis 
(CMS)  [15].  CMS  is  based  on  the  idea  that  structural  systems  can  be  modeled  using 
substructuring  techniques.  That  is,  a structure  is  comprised  of  multiple  components  and 
each  component  has  a basis  set  which  is  chosen  specifically  for  that  component. 

A most  intuitive  choice  for  the  basis  vectors  is  the  eigenvectors  of  the  component. 
Early  work  by  Hurty  [32]  and  Craig  [16]  augmented  the  modal  basis  with  additional 
modes  that  represented  the  interfacing  properties  each  component  had  with  its  inboard 
and  outboard  counterparts.  The  most  common  mode  sets  used  in  dynamic  simulations 
are  the  Craig-Bampton  modes  [16]  and  the  MacNeal-Rubin  modes  [47,  55]. 

The  Craig-Bampton  modes  are  composed  of  a set  of  fixed  interface  normal  or 
eigenmodes  and  a set  of  constraint  modes.  The  constraint  modes  are  created  by  statically 
imposing  a unit  displacement  on  one  of  the  interface  DOFs  while  keeping  the  remaining 
interface  DOFs  fixed.  The  MacNeal-Rubin  modes  consist  of  a set  of  free  interface 
normal  modes  and  a set  of  residual  attachment  modes.  The  attachment  modes  are  created 
by  statically  imposing  a unit  force  on  one  of  the  interface  DOFs  while  keeping  the 


remaining  DOFs  force  free. 
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Numerous  other  articles  have  been  written  that  expand  on  the  CMS  method.  A 
selection  of  articles  which  differ  in  applications  and  choice  of  basis  selection  can  be 
found  in  Refs.  6,  17,  76,  and  77. 

2.2.3  Modal  Ranking 

Once  a set  of  basis  modes  is  selected,  it  may  be  difficult  to  determine  if  all  the 
modes  are  really  necessary  for  an  accurate  simulation.  Some  of  the  modes  selected  may 
contribute  little  or  not  at  all  to  the  simulation  of  the  flexible  motion.  One  of  the  first 
methods  used  to  determine  modal  importance  was  initiated  by  Skelton  et  al.  [61]  who 
made  use  of  Skelton’s  component  cost  analysis  [60].  In  this  method,  each  mode’s 
contribution  is  based  upon  relative  contributions  to  a stochastic  cost  function.  The  modes 
are  then  ranked  according  to  their  contribution  and  deleted  accordingly.  One  of  the  major 
drawbacks  to  this  method  is  that  the  ranking  procedure  is  independent  of  actuator/sensor 
placement. 

An  alternative  method  for  modal  ranking  that  is  influenced  by  actuator/sensor 
placement  is  the  internally  balanced  realization  technique  which  was  first  introduced  by 
Moore  [51].  The  balanced  realization  method  is  based  on  a linear  system  being  mapped 
into  an  equivalent  system  (balanced),  where  the  controllability  and  observability 
grammians  are  diagonal  and  equivalent.  The  diagonal  terms  are  often  referred  to  as 
Hankel  singular  values  (HSVs)  and  can  be  used  to  rank  the  importance  of  the  modes. 
The  transformation  to  the  balanced  states,  however,  can  be  costly  for  large  systems  and 
there  is  a loss  of  physical  significance  to  the  state  vector.  In  a paper  by  Gregory  [25],  it  is 
shown  that  for  structural  systems  with  light  damping  and  distinct  eigenvalues,  the 
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transformation  need  not  take  place  and  the  HSVs  can  be  approximated.  Applications  of 
the  balanced  realization  methods  for  flexible  multibody  systems  can  be  found  in 
numerous  articles.  In  work  by  Spanos  and  Tshua  [65],  the  balanced  realization  method  is 
used  in  conjunction  with  a variety  of  CMS  approaches  on  a model  of  the  Galileo 
spacecraft.  Results  showed  that  all  the  CMS  methods  considered  provided  accurate 
reduced-order  models.  In  a paper  by  Kammer  and  Triller  [35],  various  modal  ranking 
techniques  including  the  balanced  realization  method,  are  compared  in  the  selection  of 
Craig-Bampton  modes.  It  was  shown  that  the  balanced  realization  method  as  well  as  the 
Effective  Interface  Mass  approaches  all  provided  efficient  methods  for  the  modal  ranking 
of  Craig-Bampton  modes.  Research  by  Gawronski  [22]  demonstrates  the  usefulness  of 
the  balanced  realization  method  in  actuator/sensor  placement. 

Fundamental  to  the  early  versions  of  CMS  is  the  selection  of  normal  modes  as  the 
basis  set.  These  modes  are  typically  chosen  to  optimally  match  the  dynamic  response 
within  a certain  frequency  bandwidth.  An  alternative  technique  for  selecting  a modal 
basis  is  based  on  parameter  matching.  That  is,  modes  are  selected  based  on  their  ability 
to  preserve  certain  parameters.  This  concept  was  introduced  in  a paper  by  Villemagne 
and  Skelton  [72],  for  the  purpose  of  generating  reduced  models  for  multi-input/multi- 
output  (MIMO)  systems.  In  their  paper,  four  parameters  are  defined  for  the  purpose  of 
parameter  matching:  low  frequency  moments,  high  frequency  moments,  low  frequency 
power  moments  and  high  frequency  power  moments.  These  parameters  are  associated 
with  the  low  frequency  response,  high  frequency  response,  low  frequency  power  spectral 
density  and  high  frequency  power  spectral  density,  respectively.  Their  approach  is  shown 
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to  yield  the  same  results  as  the  conventional  methods  while  being  more  computationally 
easy  to  obtain. 

It  has  been  shown  by  Su  and  Craig  [68],  that  with  a basis  set  defined  by  Krylov 
vectors,  particular  low  frequency  moments  are  preserved.  The  Krylov  vectors  have  a 
number  of  advantages  over  the  normal  modes.  First,  for  large  systems  they  are 
computationally  less  expensive  to  calculate.  No  eigensolution  is  required  to  form  the 
Krylov  vectors.  Secondly,  by  their  definition,  the  Krylov  vectors  contain  influences  from 
actuators/sensors.  One  drawback  to  the  Krylov  vectors  is  associated  with  their 
formulation.  Numerical  difficulties  can  arise  associated  with  their  lack  of  orthogonality 
when  a large  number  of  Krylov  vectors  are  required.  For  a complete  discussion  on  the 
Krylov  space  formulation,  refer  to  Ref.  24.  Additional  articles  demonstrating  the 
usefulness  of  Krylov  vectors  in  the  reduction  of  flexible  structures  can  be  found  in  Refs. 
17  and  67. 

2.2.4  System  Level  Component  Modes 

In  the  previous  sections,  the  selection  of  component  modal  bases  is  determined  by 
evaluating  criteria  at  the  component  level.  An  alternative  approach  is  to  consider 
evaluating  the  selection  criteria  on  modes  generated  at  the  system  level.  This  idea  was 
first  introduced  by  Eke  and  Man  [20].  At  the  system  level,  a set  of  desired  system  normal 
modes  is  projected  onto  each  component,  yielding  a reduced  set  of  component  modes  for 
each  component.  Upon  assembly  of  the  reduced  components,  the  reduced  system 
contains  the  desired  normal  modes.  Bernard  [10]  also  investigates  the  projection  and 
assembly  (P&A)  method  and  shows  why  the  desired  modes  are  retained  and  determines 
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necessary  conditions  for  a successful  reduction.  Lee  and  Tsuha  [45]  showed  that  the 
system  level  static  gain  is  not  preserved  when  only  normal  modes  are  used  in  the 
projection.  They  introduce  an  enhanced  projection  and  assembly  (EP&A)  method  where 
static  correction  modes  are  introduced  either  at  the  system  level  [45]  or  at  the  component 
level  [43]  and  the  preservation  of  the  static  gain  is  guaranteed. 

A drawback  to  the  projection  and  assembly  methods  is  that  the  reduced  system 
contains  extraneous  modes  which  lie  outside  the  bandwidth  of  the  desired  system  modes. 
In  Ref.  42,  an  extension  to  the  P&A  method  is  provided  which  identifies  and  discards  the 
component  modes  that  contribute  to  the  generation  of  the  extraneous  modes.  Additional 
work  by  Lee  and  Tsuha  [44]  determines  desired  system  level  modes  based  on  various 
configurations  of  articulating  flexible  systems. 

2.2.5  Geometric  Nonlinear  Formulations 

One  of  the  benefits  of  using  the  floating  reference  method  in  dynamic  simulation 
is  that  it  is  easily  incorporated  into  current  simulation  codes.  The  assumed  modes 
approximation  can  then  be  utilized  for  modeling  flexible  elements.  However,  in  systems 
that  exhibit  high-speed  rotations,  nonlinear  stiffness  effects  must  be  accounted  for.  Even 
for  components  operating  within  the  linearly  elastic  range,  fundamental  assumed  modes 
in  the  floating  reference  frame  cannot  model  spin  stiffening  effects.  This  modeling  error 
is  attributed  to  a premature  linearization  of  the  flexible  component’s  elastic  equations  [7]. 
Numerous  formulations  extend  the  floating  reference  method  to  capture  the  nonlinear 
stiffness  effects  [30,  36,  48,  69,  70].  In  essence,  additional  stiffness  terms  are  added  back 
to  the  component’s  linearized  equations  of  motion.  This  approach  has  two  benefits. 
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First,  since  the  local  reference  frame  is  utilized,  the  linear  FEMs  of  the  flexible 
components  can  still  be  used.  Secondly,  incorporation  of  the  nonlinear  stiffness  matrices 
into  dynamic  simulation  code  is  feasible  [30]. 

An  alternative  approach  is  the  geometrically  exact  formulation  in  which  elastic 
deformations  and  overall  gross  motions  are  referenced  to  an  inertial  frame  [62,  75].  This 
approach  employs  nonlinear  finite  element  methods  and  even  though  it  captures  nonlinear 
stiffness  effects,  it  generally  results  in  models  with  a large  number  of  DOFs.  Since  the 
formulation  is  generated  in  an  inertial  frame,  it  is  also  not  feasible  to  incorporate  the 
method  into  existing  dynamic  simulation  codes.  A thorough  review  of  the  many 
geometric  stiffening  formulations  can  be  obtained  in  Ref.  59. 

2.3  Constraint  Stabilization 

Nonlinear  mulitbody  dynamic  formulations  can  be  divided  into  two  classes, 
recursive  and  non-recursive  methods.  Recursive  formulations  generate  equations  of 
motion  using  a minimal  set  of  generalized  coordinates  which  are  defined  using  relative 
DOFs  [4,  28,  63].  Equations  of  motion  are  generated  using  a recursive  multi-pass 
procedure  in  which  kinematic  relations  are  first  defined  followed  by  the  derivation  of 
dynamical  equations.  Applications  for  recursive  formulations  are  best  suited  for  open  or 
tree  topological  type  systems  and  have  been  used  in  the  dynamic  simulation  package, 
TREETOPS  [64].  Although  recursive  methods  have  an  order  N computation  cost,  for 
rigid  body  systems,  they  present  difficulties  associated  with  implementing  constraints  and 


generalized  forces  [57]. 
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The  alternative  approach  to  recursive  methods  is  the  non-recursive  approach. 
Within  this  set  are  redundant  formulations  of  nonlinear  multibody  systems  that  generate 
dynamical  equations  in  the  form  of  DAEs.  Numerical  solutions  to  DAEs  representing 
dynamical  systems  present  a number  of  theoretical  difficulties.  In  one  class  of  solution 
techniques,  the  second  time  derivative  of  the  constraints  are  imposed  instead  of  the 
constraints  themselves.  Solutions  of  these  equations  can  be  obtained  by  numerical 
integration  in  two  ways:  (1)  by  directly  integrating  a full  set  of  equations  that  defines  the 
states  as  both  the  independent  and  dependent  coordinates  and  the  Lagrange  multipliers  or 
(2)  by  integrating  a minimal  set  or  independent  set  of  coordinates  and  then  use  the 
kinematic  constraints  to  determine  the  dependent  coordinates.  Some  of  the  latter 
formulations  include  the  range  space  and  null  space  formulations  [1,  21,  40,  46]. 
Whether  the  full  or  minimal  set  formulations  are  used  in  the  numerical  integration,  both 
methods  are  subject  to  numerical  errors  since  both  methods  impose  the  second  derivative 
of  the  constraint  equations.  These  numerical  errors  have  been  termed  constraint  violation 
drift  in  the  literature. 

There  are  a number  of  methodologies  that  seek  to  minimize  or  bypass  the  drift 
errors  that  occur  in  these  redundant  formulations.  Baumgarte  [8]  pioneered  one 
methodology  that  sought  to  stabilize  the  violation  drift.  In  his  approach,  the  constraints 
that  exhibit  the  drift  violations  are  replaced  with  a modified  set  of  constraints  that 
demonstrate  improved  numerical  accuracy.  Additional  work  has  extended  Baumgarte’s 
method  by  using  a variety  of  different  modified  set  of  constraints  [2,  5,  33].  In  work  by 
Kurdila  et  al.  [38],  a comparative  study  of  the  various  formulations  and  stabilization 
methods  are  provided  for  dynamical  systems  with  multiple  singular  configurations.  In 
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this  paper,  the  variable  sliding  mode  approach  is  shown  to  provide  the  best  stabilization 
results. 

A drawback  to  the  Baumgarte’s  stabilization  method  is  that  there  is  no  known 
procedure  that  determines  the  optimal  values  of  the  needed  parameters.  An  alternative 
method  that  provides  a more  robust  stabilization  is  provided  in  Refs.  19  and  26.  In  this 
method,  instability  is  avoided  by  the  repeated  projection  of  the  numerical  solution  onto 
the  solution  manifold. 

A class  of  solution  techniques  that  does  not  impose  the  second  derivative  of  the 
constraint  equations  attempts  to  systematically  eliminate  the  multipliers,  thereby, 
eliminating  constraint  violations  all  together.  In  work  by  Wehage  and  Haug  [27,  73],  a 
coordinate  partitioning  algorithm  is  derived  to  develop  a minimal  basis  from  partitioning 
the  Jacobian  matrix.  Although  the  minimal  set  formulations  may  appear  to  be  more 
attractive  computationally,  they  are  in  fact,  more  costly  for  some  large  systems.  This  is 
primarily  due  to  the  fully  populated  inertia  matrix  obtained  from  the  reduced  formulation 
as  well  as  the  number  of  matrix  inversions  required  to  obtain  the  reduced  formulation. 
The  dynamic  simulation  code,  DADS,  employs  a hybrid  stabilization  method  involving 
both  the  coordinate  partitioning  method  and  Baumgarte’s  constraint  stabilization  [27]. 

A final  method,  which  does  not  seek  to  eliminate  or  stabilize  constraints,  but 
approximates  the  constraints  in  a numerically  stable  manner,  is  the  penalty  method  [9, 
39].  In  this  method  the  constraint  multiplier  is  replaced  with  and  approximated  multiplier 
using  a penalty  parameter.  The  result  leads  to  an  approximate  Lagrangian  that  contains 
penalty  terms  associated  with  the  inertia,  damping  or  stiffness  of  the  system.  In  Ref.  41 
sufficient  conditions  are  derived  that  guarantee  convergence  and  Lyapunov  stability  for  a 
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class  of  multibody  systems.  Although  the  penalty  method  is  relatively  simple  to 
implement,  a major  difficulty  is  to  choose  parameters  large  enough  to  be  effective,  but 
not  so  large  to  cause  numerical  problems  [53]. 

2.4  Switched  Systems 

The  industrial  need  for  accurate  model  reduction  techniques  of  multiflexible  body 
systems  has  been  well  established  and  numerous  articles  have  addressed  these  model 
reduction  needs.  As  mentioned  previously,  one  of  the  most  popular  techniques  is  to 
model  linear  flexible  structures  using  the  assumed  modes  approach.  This  approach 
allows  for  relatively  accurate  representations  of  flexible  structures  using  a modal  basis  as 
a parameterization.  Additional  articles  use  various  modal  bases  in  • search  of  more 
accurate  responses.  Craig  augments  normal  modes  with  attachment  and  constraint 
modes.  Other  work  utilizes  Krylov  modes  as  a basis  instead  of  normal  modes.  These 
approaches  can  be  optimized  for  components  whose  boundary  conditions  and  forces  do 
not  vary  with  time.  It  is  not  the  case,  however,  for  components  that  are  subject  to  these 
conditions.  In  this  work,  switched  system  modeling  is  used  to  optimize  the  modal 
selection  procedures  as  a function  of  time. 

Switched  systems  are  systems  that  contain  both  continuous  and  discrete  elements. 
An  example  switched  system  is  one  that  comprises  a set  of  ODEs  that  switch  according  to 
a prescribed  switching  sequence.  A vast  amount  of  literature  has  been  published  that 
addresses  the  modeling  and  control  of  switched  or  hybrid  systems  [52,  56,  71,  74]. 
Stability  issues  have  also  provided  a basis  for  numerous  research  articles  to  focus  on.  In 
Refs.  12,  18,  and  54,  multiple  Lyapunov  functions  are  introduced  as  tools  for  analyzing 
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stability  in  a Lypunov  sense.  Applications  for  switched  system  modeling  can  be  found  in 
a variety  of  areas  including,  robotic  systems,  motor  dynamics,  and  computer  hardware  [3, 
13,  23].  Until  recently  [49],  however,  switched  system  research  has  not  been  applied  to 
structural  modeling  to  the  author’s  knowledge. 


CHAPTER  3 

SYSTEM  LEVEL  MODEL  REDUCTION 


3.1  Introduction 

In  this  work,  the  assumed  modes  approach  is  the  chosen  method  for  the  model 
reduction  of  linear  structural  systems.  This  approach  requires  the  selection  of  a basis  set 
of  functions  or  modes.  The  choice  of  the  basis  is  essentially  the  reduction  of  the  system 
model.  That  is,  the  dimension  of  the  basis  dictates  the  size  of  the  system.  The  fewer  the 
modes  used  to  model  the  linear  elastic  behavior,  the  smaller  the  size  of  the  resulting 
governing  equations.  Clearly,  the  choice  of  modal  basis  will  effect  the  accuracy  of  the 
simulation,  but  it  also  influences  the  cost  or  computation  time  of  the  simulation.  As  the 
number  of  kept  basis  modes  is  increased,  the  accuracy  may  increase  but  at  the  expense  of 
a higher  cost. 

The  resulting  equations  of  motion  from  the  assumed  modes  approximation 
generally  do  not  have  a closed  form  solution  and  require  numerical  methods  to  obtain  an 
approximate  solution  in  time.  The  integration  step  size  used  in  the  numerical  simulation 
is  dependent  on  the  frequencies  of  the  basis  vectors.  That  is,  the  higher  the  frequencies 
retained  in  the  mode  set,  the  lower  the  integration  time  step  must  be  and  therefore,  the 
higher  the  cost  of  the  simulation.  This  defines  the  essential  question  for  model  reduction. 
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How  does  one  choose  a minimal  basis  set  with  the  lowest  possible  frequencies  and  still 
capture  the  salient  features  of  the  dynamical  system? 

Conventional  assumed  modes  methods  use  bases  or  modes  derived  at  the 
component  level.  In  this  chapter,  modal  bases  are  determined  from  the  system  level  as 
opposed  to  the  component  level,  resulting  in  component  models  which  contain  a truer 
representation  of  inboard  and  outboard  body  effects.  The  development  of  standard 
system  level  component  modes  utilizes  P&A  concepts  and  generates  system  level  modes 
based  on  normal  modes  of  vibration.  This  work  extends  on  these  concepts  by  generating 
system  level  modes  based  on  Krylov  modes.  An  additional  extension  is  made  by 
utilizing  the  balanced  realization  theory  to  determine  the  relative  importance  of  each 
mode.  A minimal  modal  set  can  then  be  chosen  based  on  these  rankings.  Numerical 
examples  are  provided  to  demonstrate  the  effectiveness  of  various  forms  of  system  level 
modes. 

3.2  Projection  Methods 
3.2.1  Projection  and  Assembly  Method 


The  P&A  method  begins  with  the  individually  modeled  system  components.  The 
system  can  be  modeled  using  any  number  of  components,  but  for  demonstration 
purposes,  only  two  components  are  used.  The  undamped  equations  of  motion  for  the  first 
component  (component  A)  is  given  by 
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where  the  subscript  p is  the  number  of  DOFs  for  component  A,xA  is  the  displacement 

vector,  Ma  and  KA  are  the  mass  and  stiffness  matrices,  respectively,  and  Gpa  and  Hj?p  are 

the  input  and  output  influence  matrices,  respectively.  The  subscript  a represents  the 
number  of  inputs  and  the  subscript  b represents  the  number  of  outputs.  Similar  equations 
can  be  written  for  the  second  component  (component  B) 
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where  q is  the  number  of  DOFs  for  component  B. 

The  system  equations  may  be  generated  by  considering  a Ritz  transformation  that 


maps  the  vectors  xA  and  into  xn  as  follows 
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where  n = p + q - i and  i is  the  number  of  constraints.  The  matrices  PAn  (a)  and  P^n  (a) 

represent  the  displacement  compatibility  relationships  at  the  components'  interface.  The 
parameter  a indicates  that  the  system  has  a variable  configuration.  Thus,  the  resulting 
equations,  shown  below,  only  describe  the  system  at  a particular  configuration 
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where 
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Similar  expressions  can  be  written  for  the  other  system  matrices.  The  modal  equivalence 


of  Eq.  (3.2.4)  can  be  written  as 
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where  E,n  is  the  modal  coordinate  and  At m contains  the  undamped  system  eigenvalues. 

In  the  P&A  method,  only  the  modes  important  to  the  system's  input-output 
mapping  are  retained.  The  displacement,  xn , can  then  be  written  as 


Xn  = [0nk  ®n,] 
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where  &nk  and  Ont  are  the  kept  and  discarded  modes,  respectively,  and  E,k  and  are  the 
corresponding  generalized  coordinates.  The  selection  of  which  modes  to  be  retained  is  an 
important  issue  and  will  be  addressed  in  Section  3.2.4.  Note  that  since  the  system  has  a 
variable  configuration,  the  kept  modes  will  need  to  be  determined  by  analyzing  their 
importance  at  all  configurations. 

Reduced  component  models  are  generated  by  making  the  following  projections, 
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where  £ k and  %Bk  the  reduced  generalized  coordinates  of  components  A and  B, 
respectively.  Substitutions  of  Eq.(3.2.8)  into  Eqs.  (3.2.1)  and  (3.2.2)  yields  the  equations 
of  motion  for  the  reduced  component  models. 


A 
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\TJ  B 
Tqk 


KKti  + KTKKtl  = Ktg 


qk  qq  qk 


qk 


Bub 

qa  a 


(3.2.9) 
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The  reduced  components  can  be  synthesized  to  generate  a reduced-order  model  of  the 
system.  It  has  been  proven  that  the  modes,  &nk , and  the  corresponding  eigenvalues  of  the 
full-order  model  are  preserved  exactly  in  the  reduced-order  model  [10]. 

3.2.2  Enhanced  Projection  and  Assembly  Method 

Although  the  P&A  method  has  shown  the  capability  of  preserving  the  modal 
characteristics  in  reduced  models,  it  does  not  preserve  the  static  gain  of  the  system.  The 
static  gain  of  a full-order  system  may  be  defined  as 

E..  = + 0.XK)Gm  (3.2.10) 

where  Ema  is  the  m x a static  gain  matrix  and  A kk  and  A „ are  the  kept  and  discarded 
eigenvalues,  respectively.  It  should  be  noted  that  if  the  system  contains  rigid  body 
modes,  the  inverses  of  the  eigenvalue  matrices  will  not  exist  and  pseudostatic  gains  must 
be  developed  for  the  evaluation. 

The  static  gain  of  a system  is  directly  related  to  the  poles  and  zeros  of  the  system's 
transfer  function.  A change  in  the  static  gain  will  alter  the  pole  and  zero  locations  in  the 
transfer  function.  Thus,  if  the  static  gain  of  the  full-order  system  is  not  preserved  in  the 
reduced  model,  some  or  all  of  the  system's  poles  and  zeros  will  be  altered  [37]. 

The  EP&A  method  addresses  the  static  gain  issue  by  using  ideas  from  the 
component  mode  synthesis  to  augment  the  retained  mode  set  with  static  correction 
modes.  The  augmentation  can  be  done  two  ways  as  shown  in  Fig.  3.1.  It  can  be  done  at 
the  system  level  where  the  kept  modes  are  augmented  before  they  are  projected  onto  the 
components  or  it  can  be  done  at  the  component  level  where  correction  modes  are  added  to 
each  of  the  reduced  component  mode  sets. 
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Component  Level  Augmentation 


System  Level  Augmentation 


Full  Order 
System  Model 


Full  Order 
System  Model 


Figure  3.1  EP&A  Model  Reduction  Approaches 
Static  correction  modes  are  often  defined  by  attachment  modes.  There  are  three 
forms  of  attachment  modes  that  can  be  utilized:  1)  the  standard  attachment  modes,  2)  the 
inertia  relief  attachment  modes  which  are  more  consistent  with  determining  the  dynamic 
response  via  the  mode-acceleration  method  and  3)  the  residual  inertia  relief  attachment 
modes  which  combine  with  component  normal  modes  to  form  a linearly  independent  set 
of  modes  [15]. 

No  matter  which  type  of  attachment  modes  are  used,  the  choice  of  which  DOF  the 
attachment  mode  (or  modes)  is  defined  from  must  be  addressed.  For  an  n x n system  or 
component,  there  are  n choices  (DOFs)  to  define  an  attachment  mode  from.  The  choice 
of  DOF  will  depend  on  how  well  the  attachment  mode  defined  by  that  DOF  preserves  the 
m x a static  gain  matrix.  Determining  which  of  the  n reduced-order  static  gain  matrices 
are  closest  to  the  full-order  static  gain  matrix  is  not  a trivial  task.  One  approach  to  this 
problem  is  to  use  MIMO  Bode  magnitude  plots  [66].  The  underlying  principle  for 
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developing  these  plots  is  that  the  magnitude  of  the  system  transfer  function  at  any 
frequency  is  bounded  above  by  its  maximum  singular  value  and  below  by  its  minimum 
non-zero  singular  value.  Thus,  only  these  two  singular  values  need  to  be  accounted  for  in 
the  reduction  process.  The  choice  of  attachment  modes  can  then  be  based  on  which  of  the 
n reduced  MIMO  Bode  magnitude  plots  (maximum  and  minimum  singular  values) 
compares  the  best  with  the  full  MIMO  plots.  This  measure  of  comparison  is  less  rigorous 
then  determining  which  of  n attachment  modes  best  preserves  an  m x a matrix. 

For  a system  level  augmentation,  the  kept  modes,  &nk,  of  Eq.  (3.2.7)  are 
augmented  with  attachment  modes,  &na , and  the  projections  and  component  equations  are 
formed  in  the  same  manner  as  Eqs.  (3.2.8)  and  (3.2.9).  At  the  component  level,  the 
projections  can  be  rewritten  from  Eq.  (3.2.8)  as 
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xa  = W8  E 
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rAPkVt 

YBqkVl 


(3.2.11) 


where  Skk  and  EBkk  are  the  eigenvector  matrices  associated  with  the  eigenvalue  problems 


of  Eq.  (3.2.9),  respectively,  and  %k  = Skk  rjk  and  = SBkkrjB . The  matrices  yApk  and  yB 


qk 


are  normalized  such  that 
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(3.2.12) 


The  mode  sets  yApk  and  yBqk  can  then  be  augmented  with  attachment  modes  as  follows, 
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(3.2.13) 
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where  rjA  and  rfa  are  the  generalized  coordinates  associated  with  the  attachment  modes 
of  components  A and  B,  respectively,  and  v = k + a . The  equations  of  motion  for  the 
reduced  components  can  be  generated  by  substituting  Eqs.  (3.2.13)  into  Eqs.  (3.2.1)  and 
(3.2.2). 


r 

r 


(3.2.14) 


These  reduced  component  equations  can  be  recombined  to  form  a reduced  set  of  system 
equations.  As  with  the  P&A  method,  the  kept  mode  set,  &nk , and  the  corresponding 
eigenvalues  are  preserved  exactly.  However  with  the  EP&A  method,  the  static  gain  of 
the  full-order  model  has  also  been  preserved  in  the  reduced-order  model. 


3.2.3  Krylov  Projection  and  Assembly  Method 

In  the  P&A  and  EP&A  methods,  the  normal  modes  of  the  system  at  a particular 
configuration  are  projected  onto  the  full-order  components.  In  the  Krylov  method,  the 
system  Krylov  modes  are  used  as  the  projection  basis  instead  of  the  normal  modes.  An 
obvious  advantage  to  this  is  that  the  full-order  solution  of  the  eigenvalue  problem, 
required  by  the  P&A  methods,  can  be  replaced  by  the  generation  of  the  full-order  Krylov 
modes. 

When  systems  are  reduced  using  normal  modes  as  a basis,  the  optimization  is 
based  on  a specific  criteria  (i.e.,  frequency  bandwith).  The  Krylov  method  optimizes  the 
reduction  based  on  preserving  a number  of  specified  parameters  called  the  low-frequency 
moments.  For  an  undamped  system  given  by 
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Mx  + Kx  = Gu 
y = Hx 


(3.2.15) 


the  Fourier  transform  of  Eq.  (3.2.15)  and  a Taylor  series  expansion  about  a>  = 0 leads  to 
the  following  output  frequency  response 


The  terms  in  the  bracket,  H(K~l  M)‘  K~'G , are  termed  the  low-frequency  moments  of  the 
system.  As  shown  in  Ref.  [68],  a system  reduction  given  by  x = Lxr  will  preserve  i low- 
frequency  moments  if  the  span  {L}  = span  { LG } where 


The  vectors  in  La  are  termed  the  Krylov  vectors. 

As  stated  previously,  the  Krylov  vectors  (space)  may  be  used  in  place  of  the 
modal  vectors  (space)  in  the  P&A  methods.  This  equates  to  replacing  the  (j)nk  with  a 
number  of  Krylov  vectors.  Attachment  modes  can  then  be  used  to  augment  the  reduced 
component  models  as  in  the  EP&A  method.  In  the  P&A  methods  the  modes  and 
frequencies  corresponding  to  the  kept  normal  modes  were  shown  to  be  preserved  in  the 
reduced  system  [10].  With  the  Krylov  projection  method,  the  low-frequency  moments  of 
the  reduced  system,  not  the  frequencies  are  preserved. 


00 


(3.2.16) 


i'=0 


(3.2.17) 


3.2.4  Balanced  Realization  Theory 


The  determination  of  which  system  modes  are  kept  and  which  are  truncated  plays 
an  important  role  in  obtaining  the  most  accurate  and  efficient  component  models.  A 
standard  approach  might  be  to  consider  only  the  modes  within  a certain  frequency 
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bandwidth.  Depending  on  the  input-output  mapping,  some  of  the  lower  frequencies  may 
not  be  as  important  as  the  higher  frequencies,  or  possibly  some  of  the  higher  frequency 
modes  may  hold  little  significance.  Remembering  that  the  cost  of  a simulation  is  directly 
related  to  the  highest  retained  frequency,  great  consideration  should  be  given  to  the 
higher  frequency  modes  before  they  are  retained.  In  this  section,  the  balanced  realization 
approach  for  modal  selection  is  described. 

Consider  a dynamic  system  which  is  written  in  the  state-space  form  as 


x = Ax  + Bu 

y = Cx 


(3.2.18) 


The  system's  observability  grammian  W0  and  controllability  grammian  Wc  can  be  defined 
as  follows 


WQ  = f"  eAT,CTCeA‘dt 
0 Jo 

W = fV'  BBTeA?,dt 

c Jo 


(3.2.19) 


A system  is  said  to  be  internally  balanced  when  equal  weight  is  placed  with  observability 
and  controllability.  This  can  be  written  as 

W0=wc  = z (3.2.20) 


where 


Z-  diag[crl,(T2,...(Tn]  (3.2.21) 

The  <j  ' s have  been  termed  “second  order  modes”  or  Hankel  Singular  Values  (HSVs)  in 
the  literature. 

The  idea  behind  using  the  balanced  realization  theory  for  model  reduction  is  to 
bring  a model  to  a balanced  state  and  remove  the  least  significant  modes.  In  this  case,  the 
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least  significant  modes  are  determined  from  the  a' s with  the  least  magnitude.  These 
truncated  states  (balanced)  correspond  to  the  modes  with  the  least  controllability  and 
observability  for  the  given  B and  C. 

It  has  been  shown  that  any  model  with  the  same  form  of  Eq.  (3.2.18)  can  be 
brought  to  its  balanced  state  using  similarity  transformations  [51].  The  transformation, 
however,  can  be  costly  for  large  models  and  there  is  a loss  of  physical  significance  to  the 
state  vector.  Gregory  showed  that  for  systems  with  light  damping  and  distinct 
frequencies,  the  transformations  to  the  balanced  state  is  not  necessary  to  determine  the 
HSVs.  His  primary  result  is  that  a system  is  approximately  internally  balanced  with 


(4C<y/ ) 1 [bib!  {®?Ct  c( )] 


-2  T 


1-1/2 


(3.2.22) 


where  Q is  the  damping  ratio  of  the  ith  mode.  The  relative  importance  of  each  system 

mode  to  the  system's  input-output  mapping  can  be  obtained  using  this  result.  As  noted  in 
the  beginning  of  this  section,  for  articulated  multiflexible  body  systems,  the  HSV's  should 
be  evaluated  at  all  necessary  configurations.  This  is  because  the  significance  of  each 
mode  depends  on  the  configuration  of  the  system. 


3.3  Numerical  Examples 


3.3.1  EP&A  Example 

In  this  section  a numerical  example  is  presented  to  demonstrate  the  projection  and 
assembly  techniques  for  model  reduction.  Initially,  full-order  components  are  used  to 
generate  full-order  mass  and  stiffness  matrices.  Defining  the  input  and  output  system 
matrices  completes  the  full-order  system  model.  The  system  modes  can  then  be  ranked 
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2 Nodes 


Figure  3.2  Example  EP&A  System 

using  the  HSV's  and  obtain  the  reduced  components  via  the  projection  and  assembly 
method. 

For  demonstration  purposes,  an  example  is  chosen  that  does  not  consist  of 
articulating  bodies,  but  is  simply  a combination  of  multiflexible  bodies.  The  system  is  a 
cantilevered  beam  consisting  of  two  components  as  shown  in  Fig.  3.2.  Each  beam  is 
modeled  as  an  Euler-Bemoulli  beam  using  two  elements  and  three  nodes.  Each  node 
consists  of  three  DOFs;  two  translational  and  one  rotational  DOF.  Component  A is 
modeled  as  a cantilevered-free  beam  and  component  B is  modeled  as  a free-free  beam. 
The  state  vector  for  each  component  may  be  written  as 


xA=[xf  yf  0*  x2a  yA  0A]T 


B r~B  *,B  aB  vB  t,B  f)B  yb  mb  aB~\T 

To  “o  xt  y i C7|  x2  y2  tf2  J 


[x 


B 

0 


(3.2.23) 


Note  that  node  0 for  component  A has  been  eliminated  due  to  the  boundary  condition. 
The  interface  coordinates  of  each  component  must  satisfy  the  following  constraint 
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A. 

(3.2.24) 


The  inclusion  of  Eq.  (3.2.24)  with  Eq.  (3.2.23)  yields  a system  with  12  DOFs.  The 
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system  state  vector  may  then  be  written  as 


0b2]T 


(3.2.25) 


For  the  examples  in  this  chapter,  the  system  parameters  are  defined  as 


1 = 0.5 

p = 1000 

where  / is  the  element  length  (ft), 


a = 0.01  / = 8.3x1 0"6 

E = 1.5x1 04 


(3.2.26) 


a is  the  cross  sectional  area  (ft2),  I is  the  moment  of 


Table  3.1  Natural  Frequencies  of  the  Full-Order  System 


Mode 

System  Frequency  (Hz) 

1 

0.016 

2 

0.098 

3 

0.276 

4 

0.487 

5 

0.547 

6 

1.013 

7 

1.537 

8 

1.627 

9 

2.579 

10 

2.792 

11 

4.048 

12 

4.231 

inertial  (ft4),  p is  the  density  (sl/ft3)  and  E is  the  elastic  modulus  (lb/ft2).  The  12  full-order 
system  frequencies  are  given  in  Table  3.1.  Note  that  for  demonstration  purposes,  the 
system  parameters  are  chosen  to  give  frequencies  of  extremely  low  magnitudes. 

With  the  input  and  output  influence  matrices  defined  as, 


010000000000 

000010000000 

000000100000 


(3.2.27) 
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and 


001000000000 

000001000000 

000000010000 

000000001000 

000000000010 

000000000001 


(3.2.28) 


the  HSVs  can  be  calculated  via  Eq.  (3.2.22)  to  determine  the  relative  importance  of  each 


system  mode  to  the  input-output  mapping.  Results  of  these  calculations  are  shown  in 
Table  3.2.  Based  on  these  results,  only  the  first,  second,  third  and  fifth  mode  will  be 


retained  in  the  model  reduction  process. 


Table  3.2  Hankel  Singular  Values 


Mode 

HSV 

1 

14.42 

2 

2.37 

3 

0.52 

4 

0 

5 

0.26 

6 

0.04 

7 

0 

8 

0.03 

9 

0.01 

10 

0 

11 

0 

12 

0 

The  next  step  is  to  chose  an  attachment  mode  which  best  preserves  the  static  gain 
of  the  system.  In  this  study,  a system  level  augmentation  is  implemented  which  means 
there  are  twelve  choices  for  the  attachment  mode.  As  mentioned  previously,  the  MIMO 


32 


Bode  plots  give  an  indication  of  which  attachment  mode  best  preserves  the  static  gain. 
For  this  example,  the  system  with  or  without  an  attachment  mode  produces  similar  static 
gains  and  thus,  similar  MIMO  Bode  plots.  This  indicates  that  for  this  example,  it  is  not 
necessary  to  augment  the  system  with  an  attachment  mode.  The  MIMO  Bode  plot  for 


Figure  3.3  MIMO  Bode  Plot 


this  example  is  shown  in  Fig.  3.3.  The  upper  and  lower  curves  are  the  maximum  singular 
values  and  minimum  non-zero  singular  values  of  the  transfer  function,  respectively.  For 
any  choice  of  attachment  mode,  the  MIMO  Bode  plots  are  identical. 

The  EP&A  method  is  then  run  to  generate  the  reduced-order  components  and 
verify  that  the  important  system  frequencies  have  been  retained  in  the  reduced-order 
system.  The  reduced  component  and  system  frequencies  are  given  in  Table  3.3.  From 
the  HSVs,  it  was  previously  determined  that  modes  1,  2,  3 and  5 from  the  full-order 
system  model  would  be  the  modes  kept  in  the  reduction  process.  Comparing  the 
corresponding  frequencies  from  Table  3.1  with  the  first  four  frequencies  of  the  reduced 
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Table  3.3  Reduced  System  Frequencies 


Mode 

Component  A Reduced 
Frequencies  (Hz) 

Component  B Reduced 
Frequencies  (Hz) 

Reduced  System 
Frequencies  (Hz) 

1 

0.063 

0 

0.016* 

2 

0.395 

0.015 

0.098  * 

3 

1.335 

0.399 

0.276  * 

4 

3.874 

1.593 

0.547  * 

5 

1.098 

system  in  Table  3.3,  shows  that  the  frequencies  have  in  fact  been  preserved  exactly.  The 
highest  frequency  mode  is  an  extraneous  mode  and  results  from  the  reduction  process. 
The  reduced  component  models  are  now  ready  to  be  included  in  a dynamic  simulation 
package  for  an  efficient  and  effective  simulation. 

3.3.2  Krylov  Example 

For  the  Krylov  demonstration,  the  cantilevered-cantilevered  beam  shown  in  Fig. 
3.4  is  used.  The  beam  is  again  modeled  with  two  components  - each  with  three  nodes. 
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Figure  3.4  Krylov  Example  System 
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Considering  the  boundary  conditions  and  the  interface  constraint,  the  resulting  state 
vector  may  be  written  as 

[xf  yf  GA  x2a  y2A  6A  x?  yf  0ff  (3.2.29) 

The  system  parameters  are  again  defined  by  Eq.  (3.2.26)  and  the  full-order  system 
frequencies  are  given  in  Table  3.4. 


Table  3.4  Natural  Frequencies 


Mode 

System  Frequency  (Hz) 

1 

0.099 

2 

0.276 

3 

0.548 

4 

0.993 

5 

1.037 

6 

1.715 

7 

2.135 

8 

2.764 

9 

3.469 

The  generation  of  the  Krylov  vectors  requires  system  input  and  output  influence 
matrices.  For  demonstration  purposes,  a single  input  source  is  used  and  is  applied  at  the 
three  DOFs  associated  with  node  2 on  component  A.  The  output  is  simply  defined  to  be 
the  transverse  deflection  of  node  2 on  component  A.  System  frequencies  and  Bode  plots 
are  generated  for  examples  using  3,  4 and  5 Krylov  vectors.  Results  of  the  frequency 
analysis  are  shown  in  Table  3.5.  Columns  1,  3,  and  5 correspond  to  the  frequencies  of  the 
Krylov  vectors  and  columns  2, 4,  and  6 correspond  to  the  reduced  system  frequencies. 
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Table  3.5  Krylov  and  Reduced  System  Frequencies 


Krylov  (Hz) 

Red  Sys  (Hz) 

Krylov  (Hz) 

Red  Sys  (Hz) 

Krylov  (Hz) 

Red  Sys  (Hz) 

3 Modes 

4 Modes 

5 Modes 

0.099 

0.099 

0.099 

0.099 

0.099 

0.099 

0.276 

0.276 

0.276 

0.276 

0.276 

0.276 

1.315 

1.315 

0.767 

0.692 

0.548 

0.548 

1.671 

1.170 

1.037 

1.012 

2.340 

2.760 

1.037 

1.717 

2.766 

Even  though  the  frequency  analysis  for  the  Krylov  vectors  may  provide 
interesting  information,  the  more  important  information  lies  with  the  low  frequency 
moments.  Table  3.6  shows  the  low  frequency  moments  for  the  reduced  and  full-order 
model  for  each  of  the  three  runs.  It  is  easily  seen  that  for  each  case,  the  low  frequency 
moments  of  the  full-order  system  have  been  preserved  in  the  reduced  system.  As 
mentioned  previously,  it  has  yet  to  be  proven  for  the  general  case,  that  the  low  frequency 
moments  are  preserved  in  the  reduced  system  when  using  this  Krylov  model  reduction 
method. 


Table  3.6  Low  Frequency  Moments 


Full 

System 

Reduced  System 

5-Modes 

3-Modes 

4-Modes 

5 -Modes 

0.3347 

0.3347 

0.3347 

0.3347 

0.832 

0.832 

0.832 

0.832 

2.1286 

2.1286 

2.1286 

2.1286 

5.4502 

5.4502 

5.4502 

13.956 

13.956 
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The  importance  of  each  Krylov  mode  can  be  seen  be  comparing  the  frequency 
response  plots  of  the  full-order  model  with  the  various  reduced-order  models.  Figure  3.5 
compares  the  full-order  frequency  response  with  the  reduced-order  response  using  three 
Krylov  vectors.  Figure  3.6  compares  the  full  response  with  the  reduced  response  using 
four  Krylov  vectors  and  Fig.  3.7  compares  the  full  response  with  the  reduced  response 
using  five  Krylov  vectors.  It  is  easily  seen  that  a good  response  will  require  at  least  five 
Krylov  vectors  in  the  reduction  process  for  this  particular  input/output  mapping.  It 


Figure  3.5  Full  vs.  Reduced  System  - 3 Modes 


Figure  3.6  Full  vs.  Reduced  System  - 4 Modes 
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Figure  3.7  Full  vs.  Reduced  System  - 5 Modes 


should  be  noted  that  for  this  system,  only  six  linearly  independent  Krylov  vectors  may  be 
obtained. 

For  comparison  of  the  various  methods,  the  cantilevered  beam  system  is  also 
reduced  using  the  P&A  and  EP&A  methods.  The  P&A  results  are  obtained  using  four 
normal  modes  and  the  EP&A  results  are  obtained  using  four  normal  modes  and  one  static 
correction  mode.  The  frequency  responses  of  these  reductions  are  shown  in  Figs.  3.8  and 
3.9.  It  can  easily  be  seen  that  the  P&A  method  does  not  compare  well,  since  the  static 
gain  in  the  full-order  system  is  not  preserved.  The  EP&A  method,  however,  does  yield 
good  comparisons  between  the  full  and  reduced  models.  It  should  also  be  noted  that  the 
successful  Krylov  and  EP&A  runs  each  required  a total  of  five  modes.  A benefit  to  the 
Krylov  method  is  that  a full-order  eigenvalue  problem  does  not  need  to  be  solved.  For 
the  EP&A  method,  a solution  to  the  full-order  eigenvalue  problem  is  required. 
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Figure  3.8  Full  vs.  Reduced  System  - P&A 


Figure  3.9  Full  vs.  Reduced  System  - EP&A 


3.4  Summary 

The  concepts  and  derivations  of  P&A  model  reduction  methods  were  presented  as 
an  alternative  to  conventional  basis  selection  methods  for  the  assumed  modes  model 
reduction  procedure.  A variation  of  the  P&A  method,  the  EP&A  method,  was  shown  to 
extend  on  the  P&A  concepts  by  augmenting  the  mode  basis  set  and  improving  the  system 
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response  by  preserving  the  static  gain.  A new  approach  utilized  a Krylov  basis  in  lieu  of 
the  normal  modes.  One  advantage  to  this  approach  was  that  solution  of  the  full-order 
eigenvalue  problem  could  be  replaced  by  obtaining  full-order  Krylov  modes. 
Additionally,  the  balanced  realization  theory  was  incorporated  into  the  model  generation 
procedure  to  provide  a useful  method  to  rank  the  relative  importance  of  each  mode  in  the 
basis  set.  Numerical  examples  demonstrated  the  ability  for  the  P&A  and  EP&A  methods 
to  match  specified  frequencies  while  the  Krylov  method  matched  specified  low  frequency 
moments.  Frequency  response  plots  also  demonstrated  the  advantage  of  the  EP&A 
method  over  the  P&A  method  with  the  preservation  of  the  static  gain.  Similar  plots  also 
showed  that  Krylov  modes  could  be  used  to  obtain  the  same  frequency  response  as  the 
EP&A  method  without  requiring  a solution  to  the  full-order  eigenvalue  problem. 


CHAPTER  4 

CONSTRAINT  STABILIZATION 


4.1  Introduction 

A general  class  of  solution  techniques  for  multibody  systems  requires  that  the 
governing  equations  be  formulated  as  DAEs.  Solution  of  the  DAEs  presents  a number  of 
theoretical  difficulties.  In  dynamical  systems,  the  difficulties  often  lead  to  the 
phenomena  known  as  constraint  violation  drift.  Numerical  methods  to  attenuate  the 
effects  of  constraint  violation  drift  have  come  to  be  known  as  constraint  stabilization 
methods  in  the  literature.  Various  authors  have  proposed  specific  methods  to  stabilize 
constraint  violation  that  make  explicit  use  of  the  special  structure  of  the  governing 
equations 

In  this  chapter,  it  is  shown  that  a nonlinear  control  theoretic  framework  provides  a 
useful  way  of  interpreting  some  of  the  conventional  constraint  stabilization  methods. 
First  it  is  shown  using  Lie  algebraic  methods,  that  a class  of  stabilization  techniques  can 
be  recast  into  a form  that  can  be  viewed  as  an  asymptotic  output  stabilization. 
Additionally,  it  is  shown  that  the  method  of  conservation  laws  for  nonlinear  control  can 
be  used  to  represent  some  penalty-based  constraint  stabilization  methods.  This  leads  to  a 
framework  that  enables  the  introduction  of  new  variations  of  constraint  stabilization.  The 
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chapter  is  concluded  with  numerical  examples  of  the  stabilization  methods  and  a 
summary  of  the  chapter. 

4.2  Dynamic  Formulation 


In  this  dissertation,  only  nonlinear  multibody  systems  who's  governing  equations 
can  be  written  as 


d 

dt 


dT\ 


dq 


dT  dV 
+ 


K^k  J 


dqk  d(ik  «•= i d(i 


D r)(k 

+ Z -^LXi=Qk  k = \,...,N 


(4.2.1) 


are  considered.  In  these  equations  T is  the  kinetic  energy,  V is  the  potential  energy, 
q{---qN  are  the  generalized  coordinates,  QX---QN  are  the  generalized  forces  and 

/l,  ••  -Ad  are  the  Lagrange  multipliers.  The  multibody  system  is  subject  to  D holonomic 


constraints  that  can  be  written  as 


0i(q(t))  - 0 i = l--  D (4.2.2) 

Differentiating  Eq.  (4.2.2)  twice  with  respect  to  time  and  coupling  it  with  Eq.  (4.2.1) 
leads  to 


[M(<7(0)]  [B(q(t))f  |<7(0l  = ip(q,q,t)  1 
[5(^(0)]  0 JU(0j  \r(q,q,t)  J 

where  M(q)  is  an  NxN  generalized  mass  matrix, 


d<t>{q(t)) 

dq 


(4.2.3) 


(4.2.4) 


is  the  DxN  constraint  Jacobian  matrix,  p(q(t),q(t),t)  eRN  and  r(q(t),q(t),t)  eR° . The 
vector  r(q,q,t)\s  obtained  by  differentiating  the  constraints  to  obtain 


[B]q  = ~[B]q  = r(q(t),q(t),t ) 


(4.2.5) 
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To  study  the  stabilization  of  constraint  violation,  the  governing  equations  in  (4.2.1)  are 
cast  in  first  order  form.  The  2N  state  vector  is  defined  as 


(4.2.6) 


and  the  governing  equations  may  be  written  as 


7=1 


= F(X(t))  + G(X(t))A 


(4.2.7) 


where 


F(X(  0)  = 
G(X(t))  = 


-[M{x^(t))Y  p{xx,x2,t) 
0 

-[M(af1(/))]",5r(jf1)J 


1 


(4.2.8) 


Numerical  methods  for  the  solution  of  the  coupled  DAEs  appearing  in  Eqs.  (4.2.2) 
and  (4.2.3)  have  been  studied  by  many  authors,  (see  Sec.  2.5  for  the  literature  review)  In 
this  dissertation,  the  problem  is  recast  to  find  a solution  X(t)  and  control  u{t)  that  satisfy 


X(t)  = F(X(t))  + G(X(t))u(t)  (4.2.9) 


It  should  be  noted  that  in  general,  the  control  u(t)  in  Eq.  (4.2.9)  does  not  equal  the 
Lagrange  multiplier  appearing  in  Eq.  (4.2.7).  The  performance  of  a specific  control  u 
will  be  measured  in  terms  of  the  following  observations 


!>,(*,(  0)1 


Y(t ) = H(X(t))  = 


l(0). 


(4.2.10) 
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An  approximate  Lagrange  multiplier  is  sought,  u(t)  ~ A(t) , such  that  = 0 . 
Even  though  the  solution  is  not  exact,  the  stability  and  robustness  properties  associated 
with  the  violation  drift  are  hopefully  improved. 

4.3  Lie  Algebraic  Control  Methods 


Lie  algebraic  methods  have  become  a vital  tool  in  the  development  of  certain 
classes  of  nonlinear  control  strategies.  These  nonlinear  methods  require  that  the  original 
problem  be  rewritten  using  the  Lie  algebraic  approach.  In  this  section,  several  definitions 
are  introduced  that  enable  the  multibody  equations  to  be  interpreted  as  a Lie  algebraic 
control  problem. 

For  any  function  V:RN  —>  R , the  Lie  derivative  of  V along  f:RN—>R  is  defined 

to  be 


JL  flV 

' ttSx, 


(4.3.1) 


Eq.  (4.3.1)  can  easily  be  thought  of  as  a directional  derivative.  First,  consider  the 
following  system  as  a single-input/single-output  (SISO)  system 


x(t)  = f(x(t ))  + g(x)u(t) 
y(t)  = h(x(t )) 


where  f:R->R,  g:R^>  R and  h.R^>  R.  The  SISO  system  has  relative  degree  r at  x0 
provided  that 

(i)  LgLkf  h(x0)  = 0 for  all  x enbhd(x0)  VA:<r-l 

(ii)  LgLrj'h(x0)  * 0 
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Alternatively,  the  relative  degree  can  be  viewed  as  "the  minimum"  number  of  times  the 
output  y(t)  must  be  differentiated  before  the  input  u{t)  appears  explicitly.  For  example, 
with  the  output  defined  as 

j>(0  = h(x(t))  (4.3.3) 

the  first  derivative  of  the  output  is 


y(0  = Z 

i=l 

N 

-z 


Sh(x(t)) 

dXj 


r)h 

— ( fi  WO)  + gi  WOMO) 

oxi 


= Lfh(x(t))  + Lgh(x(t))u(t) 


= Lfh(x(t )) 


Similarly,  the  additional  derivatives  are 


(4.3.4) 


(4.3.5) 


/r)(0  = L{ph(x(t))  + LgL(rf-x)h{x{t))u{t) 

A standard  result  of  Lie  algebraic  control  is  that  if  the  relative  degree  r = N,  then  it  is 
possible  to  completely  feedback  linearize  the  SISO  system.  For  this  case,  it  is  possible  to 

define  a new  set  of  variables  as 


€i(t)  = h(x(t)) 

£2(0  = Lfh(x(t)) 

. ' (4.3.6) 

ZAO  = L(/-1)  h(x(t)) 


By  construction,  the  control 
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u(t)  = [LgLr-'h{x{t))Y  (-Uph(x(t))  + v(0) 


(4.3.7) 


is  feedback  linearizing.  This  means  if  u(t)  in  Eq.  (4.3.7)  is  substituted  into  (4.3.2)  a linear 
system  with  a new  control  input  v{t)  is  obtained. 

The  above  methodology  can  be  extended  to  MIMO  systems  where  the  number  of 


inputs,  D,  is  equal  to  the  number  of  measurements  or  constraints.  Differentiating  the  p 


th 


output  until  the  input  u appears  leads  to  the  following  general  expression: 


Yp(t)  = Hp(x(t)) 

JL,  dH  (x(t)) 

y,w  = E — — mo 


i= 1 


dx :. 


N dH 

= ^^{F(x(t))  + G(x(t))u(t)] 

1=1  VX- 


(4.3.8) 


D 


= LFHp  (x(t))  + £ LCj  Hp  (. x(t))Uj  (t) 


7=1 


= LFH  p{x(t)) 


7=1 


= L2FHp(x(t )) 


(4.3.9) 


y/'-’w  = L$%(m)f.Lo,(0~UH,)uAO 

7=1 

As  in  the  case  of  the  SISO  system,  a new  set  of  coordinates  corresponding  to  the  p,h 
output  are  defined  as  follows 


' ' 

> 

-J 

< . 
• 
• 

> — < 

lfhp 

• 

5 P'rr  . 

r 

1 

v 

(4.3.10) 
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Differentiating  the  new  coordinates  allows  for  a new  set  of  governing  equations  to  be 
written  as 


f,.l 

%p,l 

• 

• 

• 

V 

II 

^ A 

j’W  . 

k 

£p,  2 


7=1 


(4.3.11) 


where  r is  the  p‘h  relative  degree.  Finally,  the  last  row  from  each  block  of  coordinates, 
corresponding  to  a particular  output  is  extracted  to  obtain  a nonlinear  system 


£l,r, 

r 

Vi 

• 

< : 

• 

pD'r° . 

> — < 

• 

• 

• 

7(ro)  H 

L^F  11  d 

V*  J 

> + 

La,  Lf 


(o-i) 


H. 


Lr.  L(:d~x)H 


Gx  F 


D 


LclV'H, 


LcA°~"Ho 


U] 


u 


D 


V J 


(4.3.12) 


The  relationship  between  feedback  linearizing  controllers  and  some  constraint 
stabilization  methods  can  now  be  made  precise.  If  the  governing  equations  have  the  form 
in  Eq.  (4.2.7)  and  (4.2.10),  differentiating  Eq.  (4.2.10)  gives 


Yp{o  = £^,(0 


w dxi 


dHp  dHp 


-if 


Y(t)  = 


dx]  dx 

8H 


2 

-1 T 


{E(x(0)  + G(jc(0M0} 


dx. 


0 


-[A/(*,  (/))]"' p(x,,x2,0 
0 

BT(xx{tj) 


+ 


u(t ) 


[#(*,  (0)]*2  + X 

LaHu(l ) 


L"'H 


(4.3.13) 


Since  the  control  u{t)  does  not  appear,  a second  derivative  is  taken 
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d ([5(x, (oj*2 ) -^-([^i(oK) 


dx] 


dx. 


0 


-[M{xx )]"'  p\  +{-[M(*. )]  1 [B(xi  (0] 
= [#(*,  (0)]*2  -[5(xi(0)][M(j:,)]  1 p 


■ u(t ) 


— v 

L'Ph 


[B(x^  (t))][M(xl )]  ‘ [fl(x,  (t))]Tu(t) 


~v 

rd) 


LGLp'H  u(t) 


(4.3.14) 


A feedback  linearizing  controller  can  be  defined  if 

LGL("H  = [B][My'[B]T 


(4.3.15) 


is  invertible.  This  matrix  is  exactly  the  coefficient  matrix  in  Eq.  (4.3.12)  that  must  be 
inverted  to  calculate  a feedback  linearizing  controller 


2,  •••  /„.  L"JH 


LcA'H< 


(1) 


'G,  D 


gd  f “D 


= lgl(Ph 


(4.3.16) 


The  feedback  linearizing  controller  is  just 


u{t)  = (LGl}"H)-\-L{?H  + v) 

= ([5(x,(0)][M(x,)]  ‘[5(x,(0)]r)  (-[5(x,(0)]^2 +[5(x,(0)][M(x,)]  ‘Jp  + v(0) 


(4.3.17) 


The  choice  of  v(t)  determines  the  stabilization  method  under  consideration.  For  example, 
Baumgarte's  method  defines  v(t)  as 


v(t)  = 2 + CD2<f)(i) 


(4.3.18) 


Similarly,  Chiou's  sliding  mode  approach  defines  the  new  controller  as 


v(0  = if, \t ) + c,^(0  <f)(t)  + c2  (ifi(t)  + c,^(0) 


(4.3.19) 
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This  framework  allows  for  the  easy  construction  of  numerous  families  on  stabilization 
methods. 

4.4  Conservation  Laws  for  Constraint  Stabilization 


In  the  previous  section,  constraint  stabilization  methods  were  defined  in  the  Lie 
algebraic  control  framework.  In  this  section,  other  constraint  stabilization  methods  are 
cast  in  terms  of  nonlinear  controls  derived  from  conserved  quantities.  It  will  be  shown 
that  the  stiffness,  damping  and  inertial  penalty  methods  can  be  construed  as  constraint 
stabilization  based  on  conservation  laws. 

Rather  than  solve  the  original  system  of  DAEs  that  include  the  Lagrange 
multipliers,  another  family  is  obtained  when  penalized  versions  of  the  kinetic  and 
potential  energies  as  well  as  the  Rayleigh  dissipation  functions  are  introduced 


T = T + 


1 i>T» 


2 £ 

1 


V = V + — (j>  a<t> 


/?  = 


2s 

R + — (/>  $ 

2s 


(4.4.1) 


In  these  equations,  T is  the  kinetic  energy,  V is  the  potential  energy,  (j>  are  the  functions 
defining  the  holonomic  constraints  and  e is  a control  parameter.  The  approximating 
governing  equations  derived  from  the  penalized  Lagrange  equations  is 


dT  A 


dt 


dq 


dT  dV  1 
+ + 


V'-'v*  j 


dq.  dqk  e 
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{0<t>  + M0  + a<f\  = Qk 


(4.4.2) 


The  matrices  \/3t\  [a,>]  and  [a,y]  are  assumed  to  be  symmetric  and  positive  definite.  It 


should  also  be  noted  that  the  holonomic  constraints  are  not  enforced  exactly  in  this 
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formulation,  as  they  were  in  the  last  section.  However,  in  some  cases  it  is  possible  to 
derive  bounds  on  the  magnitude  of  constraint  violation  as  shown  in  Ref.  [41]. 

While  it  may  not  be  immediately  apparent  that  the  choice  of  nonlinear  control  in 
Eq.  (4.4.2)  yields  an  approximation  to  the  constrained  multibody  system,  it  is  true  that  as 
s — > 0 the  solution  of  Eq.  (4.4.2)  approaches  the  solution  of  the  original  DAEs,  provided 
some  common  restrictions  hold.  For  example,  if  the  system  is  natural  (T=T2), 

ju  = 0 and  (9  = 0,  it  possible  to  show  that  there  is  a Jacobi  integral  for  the  system.  This 
conservation  law  takes  the  form 


dt 


(7’+F)  + ^|-{«(’'^  + fJra(S}  = 0 


(4.4.3) 


Since  the  initial  conditions  satisfy  the  constraints  exactly,  that  is 


mm = o 

fcq{0),qm  = 0 


(4.4.4) 


integration  of  Eq.  (4.4.3)  yields 


. -r  R . 1 T (X 

T+V  + + -0t-<P  = Tq+Vo=Eo 

2 s 2 e 


(4.4.5) 


where  the  constant  E0  does  not  depend  on  s.  It  can  easily  be  seen  that 


m\f  + ||(*<)E  s se0 


(4.4.6) 


where 


Afp  = ATpA. 
A 2 = ATaA 

a 


(4.4.7) 
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for  any  vector  A.  Thus,  the  expression  in  Eq.  (4.4.6)  provides  a uniform  bound  on  the 
constraint  violation.  In  fact,  it  is  possible  to  show  that  the  solution  pair  (qe,q£)  of  Eq. 
(4.4.2)  approaches  the  solution  ( q,q ) of  the  original  system  of  DAEs. 


as  (4.4.8) 

It  should  be  noted  that  violation  bounds  for  the  Lie  algebraic  control  methods  in  the  last 
section  were  not  obtained. 

The  approach  described  above  relied  on  the  choice  of  nonlinear  control 


u = 


• • • 

P(j)  + iu(j)  + a<f> 


(4.4.9) 


in  Eq.  (4.4.2).  This  choice  enabled  the  derivation  of  the  Jacobi  integral  in  Eq.  (4.4.3). 


More  generally,  it  is  possible  to  show  that 


—(T+V)  = -'fu  (4.4.10) 

dt 

It  remains  an  open  question  whether  a choice  of  u corresponding  to  a sliding  mode 
control,  as  derived  by  Chiou  in  the  Lie  algebraic  formulation,  will  be  effective  in  this 
context. 


4.5  Numerical  Examples 


The  numerical  example  studied  in  this  section  is  depicted  in  Fig.  4.1.  The 
example  has  been  selected  to  illustrate  how  multibody  sub-assemblies  can  be  constrained 
via  the  control  theoretic  frameworks  presented  in  this  paper.  Two  sub-assemblies,  shown 
in  Fig.  4.2,  comprise  the  overall  model:  (1)  the  base  vehicle  assembly,  and  (2)  the 
weapon  platform  assembly.  The  base  vehicle  assembly  is  a full  vehicle  model  for  the 
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Figure  4.1  HMMWV/Weapon  Assembly 


highly  maneuverable  multi-wheeled  vehicle  (HMMWV)  that  has  been  developed  at 
Redstone  Arsenal  by  Dr.  Mike  Hale  and  his  research  associates.  The  base  vehicle  model 
is  comprised  of  36  rigid  bodies  represented  in  terms  of  252  generalized  coordinates. 
Because  the  orientation  of  each  body  is  expressed  in  terms  of  Euler  parameters,  there  is 
one  Euler  parameter  constraint  per  body  for  a total  of  36  Euler  parameter  constraints.  In 
addition  there  are  6 ground  constraints,  12  revolute  joins,  16  spherical  joints,  and  17 
brackets.  In  all,  there  are  consequently  212  scalar  constraints  imposed  in  the 
representation  of  the  base  vehicle  model.  When  redundant  constraints  are  eliminated,  the 
base  vehicle  has  19  DOFs.  Specific  properties  of  the  HMMWV  model  are  proprietary 
information  and  are  left  out  of  this  document. 

The  weapon  platform  subassembly  has  been  derived  in  consultation  with  Mike 
Mattice  of  Picatinny  Arsenal.  It  is  comprised  of  three  bodies  with  sufficient  constraints 


so  that  only  elevation  and  azimuth  motions  are  permitted. 
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Figure  4.2  Sub-Assemblies 

The  simulations  are  run  using  the  computer  software  package,  DADS.  The 
HMMWV  is  prescribed  to  move  forward  at  10  mphs  and  is  subjected  to  input 
disturbances  in  the  form  of  2 inch  bumps  at  intervals  of  every  2 feet.  The  approximate 
constraints  between  the  two  sub-assemblies  are  imposed  between  the  four  comers  of  the 
base  plate  of  the  weapon  assembly  and  the  top  of  the  HMMVW  cabin. 

Figure  4.3  shows  the  constraint  violation  at  one  of  the  comers  of  the  assemblies 
using  damping  and  stiffness  penalty  methods.  An  initial  value  of  500,000  lb/in  is  chosen 
for  the  stiffness  and  zero  for  the  damping.  As  the  stiffness  is  increased  by  an  order  of 
magnitude  to  5,000,000  and  then  again  to  50,000,000,  the  violation  decreases  by  an  order 
magnitude  from  0.001  to  le'4and  then  to  le'5.  With  damping  set  at  5,000  lb/in/sec,  each 
of  the  three  previous  stiffness  methods  remain  at  the  same  order  of  violation,  as  expected. 


but  have  smaller  and  smoother  oscillations. 
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Figure  4.3  Penalty  Method  - Position  Violation 


Figure  4.4  shows  the  corresponding  constraint  violation  in  derivatives  (j){t)  of  the 

simulations  described  above.  Here,  the  velocity  constraint  violation  cf>{t)  is  only  slightly 
improved  as  the  stiffness  is  increased  and  the  damping  is  set  to  zero.  With  damping 
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Figure  4.4  Penalty  Method  - Velocity  Violation 
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included,  an  improvement  of  an  order  of  magnitude  is  demonstrated  and  improves  as  the 
stiffness  increases. 

The  next  two  figures  summarize  results  obtained  with  the  hybrid  penalty  method. 
This  method  is  derived  from  the  penalty  methods  and  uses  a coefficient  defined  by  a 
sliding  mode  law.  The  hybrid  method  employs  a feedback 

w(0  = c(t)<j){t)  + K<f>(t)  (4.5.1) 

where 


m + Dm 


s 


max 


(4.5.2) 


For  the  hybrid  method,  the  stiffness  coefficient  is  set  to  500,000  lb/in  and  the 
proportionality  constant,  D,  is  set  to  1.  Three  cases  are  run  with  a Smax  that  decreases 

from  1/100,000  to  1/10,000,000.  Figure  4.5  shows  the  position  violation  <p{t ) for  each 


Figure  4.5  Hybrid  Method  - Position  Violation 
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Figure  4.6  Hybrid  Method  - Velocity  Violation 
value  of  Smax  and  Fig.  4.6  shows  the  corresponding  velocity  violations  <f>(t) . There  seems 
to  be  very  little  improvement  in  the  violations  when  Smax  is  decreased  an  order  of 
magnitude  from  1/100,000,  however,  there  is  a noticeable  improvement  when  Smax  is  set 
to  1/10,000,000. 

Figure  4.7  shows  a comparison  of  the  velocity  violations  for  both  the  penalty  and 
the  hybrid  methods.  In  Fig.  4.7  two  penalty  method  cases  are  superimposed.  The  two 
penalty  method  cases  have  the  stiffness  coefficient  set  to  be  500,000  lb/in,  which  is  the 
same  for  the  hybrid  cases,  and  the  damping  values  are  set  at  5,000  and  50,000  lb/in/sec. 
Thus,  Fig.  4.7  is  a comparison  between  using  constant  damping  coefficients  (penalty 
method)  and  variable  damping  coefficients  (hybrid  method).  The  results  indicate  that  the 
penalty  method  with  c=5000  lb/in/sec  offers  only  a slight  improvement  in  the  velocity 
violation  over  the  first  two  hybrid  methods.  Similarly,  the  penalty  method  with  c=50,000 
lb/in/sec  offers  only  a slight  improvement  over  the  most  accurate  hybrid  method. 
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Figure  4.7  Velocity  Violation 
4.6  Summary 

In  this  chapter  it  was  shown  that  the  constraint  stabilization  problem  can  be 
viewed  as  a nonlinear  controls  problem  using  two  different  formulations.  With  the  use  of 
Lie  algebraic  methods,  the  problem  can  be  cast  as  a feedback  linearizing  controller  and 
the  classic  stabilization  methods  from  Baumgarte  and  Chiou  are  easily  obtained. 
However,  if  the  problem  is  formulated  using  energy  principles,  the  penalty  methods  can 
be  obtained.  These  alternative  frameworks  lend  to  the  development  of  and  extension  to 
new  stabilization  methods.  Numerical  examples  that  define  the  stabilization  between  two 
subsystems,  a HMMWV  and  a weapon  assembly,  are  provided.  Results  indicate  that  a 
hybrid  method  described  as  a combination  of  the  penalty  methods  and  the  variable  sliding 
mode  method  provides  a competitive  alternative  to  the  classical  methods  of  stabilization. 


CHAPTER  5 
SWITCHED  SYSTEMS 


5.1  Introduction 

An  alternative  approach  to  the  conventional  methods  of  model  reduction  and  the 
methods  discussed  in  Chapter  3 models  the  structures  as  switched  systems.  In  essence, 
this  means  that  at  different  times,  the  model  of  the  system  is  changed  according  to  a 
switching  sequence.  As  boundary  forces  on  a component  vary  with  time,  the  component 
model  can  be  switched  to  contain  a basis  that  is  relevant  for  the  current  boundary  forces. 
In  this  work,  at  each  switch  time,  the  component  modes  are  ranked  according  to  modal 
participation  factors,  or  error  indicators.  The  least  important  modes  are  removed,  and 
alternative  selections  of  modes  are  inserted  in  the  active  set.  Thus,  a switched  dynamical 
system  is  obtained. 

There  are  a number  of  issues  to  address  when  modeling  structural  systems  as 
switched  systems.  How  often  the  system  is  switched,  what  logic  determines  the  number 
and  choice  of  modes,  and  the  determination  of  stability  are  just  a few  questions  which 
will  need  to  be  answered.  As  a first  step,  only  linear  structural  systems  are  considered  for 
the  switched  system  modeling.  In  this  chapter,  a general  discussion  on  switched  systems 
is  given  and  a switched  model  is  derived.  This  is  followed  by  a general  discussion  on  the 
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stability  or  instability  of  dynamical  systems.  Stability  conditions  for  the  proposed 
switched  model  are  then  derived. 

5.2  Switched  System  Model 

5.2.1  Model  Definition  and  Switching  Logic 

This  dissertation  will  consider  a switched  system  that  is  defined  as 

MO  = ft(MO)>  i eQz  (5.2.1) 

where  x{t ) € R" . Each  of  the  f’s  must  be  globally  Lipschitz  continuous  and  the  fs  must 
be  defined  in  such  a way  that  there  are  finite  switches  in  finite  time.  Additional 
constraints  can  be  added  for  systems  that  require  continuous  states.  These  systems  are 
called  continuous  switched  systems  and  have  the  added  requirement  that 

fkM  (x(tj),tj)  = fk,  ( x(tj),tj ) (5.2.2) 

The  first  task  will  be  to  cast  the  linear  structural  dynamics  equations  in  the  form  shown  in 
Eq.  (5.2.1).  Consider  a dynamical  system  of  the  form 

[m]x(t)  + [c]  x(0  + [k]x(t)  = f(t)  (5.2.3) 

where  m,  c,  and  k e Rnxn  are  the  mass,  damping,  and  stiffness  matrices,  respectively.  The 

vectors  x(t)  and  ft)  e R”x'  are  the  state  vector  and  forcing  function,  respectively.  In  state 
space  form,  the  equations  are 

Z(t)  = [A]Z(t)  + F(t)  (5.2.4) 

where  Z(t)  and  F(t)  e R2"''  are  the  new  state  and  forcing  vectors.  In  terms  of  model 
reduction,  a standard  approximation  for  the  given  system  is  the  Ritz  approximation.  The 


Ritz  approximation  takes  the  form 
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Z(t)  = M 77(0  (5.2.5) 

where  y/e  R2nxm  is  a matrix  whose  columns  are  comprised  of  the  Ritz  vectors  and 

J2(t)  € R mx 1 is  the  vector  of  generalized  coordinates. 

In  this  dissertation  an  alternative  form  of  the  Ritz  approximation  is  used; 

Zi(t)  = [iRi]rji^)  + Z0i  (5.2.6) 

where  i refers  to  the  ith  switch  and  Z0j  is  the  initial  state  vector  for  the  i‘h  model  used  after 

the  i,h  switch.  In  the  remainder  of  this  article  the  phrase  “i,h  model”  will  be  used  to  denote 
“the  model  used  after  the  i‘h  switch  time.”  Substitution  of  Eq.  (5.2.6)  into  Eq.  (5.2.4) 
generates  the  reduced  equations 

7|.(0  = [^Jr[^][^j7,.(0  + [^Jr([^M0i+Z(0)  (5-2.7) 

Eqs.  (5.2.6)  and  (5.2.7)  combine  to  form  the  set  of  dynamic  equations  for  the  i'h  model. 

rp  9 

Note  that  the  derivation  of  Eq.  (5.2.7)  has  enforced  the  conditions  [y/t]  [^,]  = / . This 
simply  means  that  the  Ritz  basis  for  the  ilh  model  must  be  orthonormal. 

The  states  Z(t ) will  be  continuous  if  the  following  conditions  are  included  as 
constraints  on  Eq.  (5.2.6). 

Z0i=Zfi]  and  77,('o)  = Q (5.2.8) 

Here,  Zf  is  the  final  value  of  the  state  vector  for  the  (i-l)lh  interval.  Evaluating  Eq. 
(5.2.6)  at  t0  and  substituting  the  conditions  in  Eq.  (5.2.8)  into  Eq.  (5.2.6)  clearly  shows 
that  state  continuity  is  preserved  between  switched  models. 

Z,(t  0)  = [V>]%(t0)  + z/M  = zfix 


(5.2.9) 
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The  basis  Yi  is  selected  as  a subset  from  a library,  IP,  which  is  determined  a 
priori;  i.e.,  c IP.  As  the  conditions  in  the  system  vary,  the  y,  s used  will  vary  within 

the  span  of  the  vectors  in  I P The  basis  library  IP  can  be  constructed  from  any  choice  of 
Ritz  vectors,  i.e.,  eigenvectors,  krylov  vectors,  constraint  vectors  or  any  combination 
thereof. 

The  most  challenging  aspect  of  this  work  is  the  design  of  a method  to  select  Yi 

from  I P Any  reasonable  method  to  select  Yi  fr°m  ^ requires  some  measure  of  modal 
participation.  As  noted  in  Section  2.2.3,  a number  of  modal  ranking  methods  have  been 
investigated.  The  choice  of  modal  ranking  in  this  work  stems  from  a portion  of  the 
reduced  governing  equations  found  in  Eq.  (5.2.7).  The  term  of  interest  is 

W i]\[A]Z0i  + Fit))  (5.2.10) 

This  term  represents  a measure  of  orthogonality  between  the  chosen  basis  for  the  i,h 
model  and  a combination  of  the  forcing  function  and  the  initial  conditions  of  the  i,h 
model.  If  Yi  is  completely  orthogonal  to  [A\Z0i  + F(t) , then  it  is  clearly  not  a good 

choice  for  a basis  in  the  i'h  interval. 

Let  Rj(t)  be  defined  as 

Ri(t)  = [A]Z0i+F(t)  (5.2.11) 

Since  Rj(t)  e Rlnx> , then  [y.YKM)  e R"‘xl  where  m is  the  number  of  modes  kept  in  y i ■ 
A scalar  quantity  representing  the  orthogonality  of  Yi  wdh  respect  to  R^t)  can  be  written 
as 


ei{t)  = Ri{t)T{Yi}[Yi}TUt) 


(5.2.12) 
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If  \f/i  spans  the  entire  space,  then  [^,][V,]  equals  the  identity  matrix,  and  the  error 
indicator  is  just  equal  to  the  norm  of  the  residual.  This  implies  that  the  maximum  value 
of  et(t)  may  be  written  as 

emaXi(t)  = Ri(t)TRi(t)  (5.2.13) 

Thus,  a measure  of  modal  participation  for  the  ith  interval  or  model  is  defined  as 

Ei  = ei  (0  / emax.  (0  (5.2.14) 

The  maximum  value  of  Ei  is  one  and  the  minimum  value  is  zero.  For  the  remainder  of 
this  paper,  Ei  will  be  referred  to  as  the  ith  error  indicator. 

With  a measure  of  modal  importance  defined,  an  algorithm  defining  the  details  of 
the  modal  selection  process  needs  to  be  developed.  Before  the  procedure  can  be 
implemented,  an  acceptable  level  needs  to  be  determined  for  the  error  indicator  in  Eq. 
(5.2.14).  Thus,  the  variable  Emin  is  introduced  and  is  defined  to  be  the  minimum 
acceptable  value  for  the  error  indicator.  It  should  be  noted  that  the  higher  Emin  is  set,  the 
more  accurate  the  simulation  will  be.  Of  course,  the  simulation  will  then  be  more 
computationally  costly.  Once  Emin  is  chosen,  the  algorithm  for  choosing  appropriate  Ritz 
vectors  must  be  developed.  To  simplify  the  algorithm,  a basis  library  E is  chosen  to  be 
a basis  for  R2n . Since  the  must  be  an  orthonormal  basis,  the  basis  library  E is 
orthonormalized.  This  means  that  for  any  y/t  e E that  is  selected,  if/j  is  guaranteed  to  be 
an  orthonormal  set.  In  general,  the  procedure  goes  as  follows: 

1 ) At  each  switch  time,  values  for  Ej  are  calculated  for  each  mode  shape 
contained  in  E That  is,  if  E e Rlnxm , then  m values  for  Ei  are  calculated. 

( 1 Ei m Ei ) The  presuperscript  indicates  which  modal  vector  is  being  used 
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from  and  the  subscript  i represents  which  model  or  switch  interval  is  being 
analyzed. 

2)  The  mode  associated  with  the  largest  value  of  Ei  is  retained  in  . 

3)  If  the  largest  value  of  Ei  is  below  Emin , then  the  mode  associated  with  the 
second  highest  value  for  Ei  is  also  retained  in  \j/i . 

4)  Number  3)  is  repeated  until  Et  > Emin . If  this  condition  cannot  be  met  with  a 
reduced  model,  then  the  entire  space  is  used  in  the  simulation  for  the  i'h 
switch  interval. 

5.2.2  Stability  Issues 

In  contrast  to  conventional  linear  time  invariant  systems,  the  stability  of  switched 
dynamical  systems  can  be  difficult  to  establish.  In  fact,  some  results  regarding  the 
stability  of  switched  systems  are  quite  counter-intuitive,  as  the  following  example 
demonstrates. 

Consider  a switched  dynamical  system  that  is  constructed  from  two  classical 
oscillators.  Let  the  switch  state  be  denoted  by  i,  and  define  a differential  system  for  i = 1 
by 

+ c,i,(t)  + k\Xx(t)  = 0 (5.2.15) 

subject  to  appropriate  initial  conditions.  Similarly,  for  the  switch  state  i = 2,  define  a 
differential  system 

m2x2(t)  + c2x2(t)  + k2x2(t)  = 0 (5.2.16) 

again  subject  to  appropriate  initial  conditions.  The  two  differential  equations  in  Eqs. 
(5.2.15)  and  (5.2.16)  can  be  used  to  define  a switched  dynamical  system  given  as 


ii„  (0  = K k it  (0 

h = *(**-1  (0) 


(5.2.17) 
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where  s is  the  switching  function.  Carefully  note  that  if  all  of  the  constants 

mi,ci,ki  g R+ , i = 1,  2,  etc.,  then  either  of  the  individual  systems  in  Eqs.  (5.2.15)  and 
(5.2.16)  is  asymptotically  stable  to  the  origin.  However,  it  is  relatively  simple  to 
construct  switching  laws  ik  = s(it_nzt_t(t))  for  the  system  in  Eq.  (5.2.17)  that  alternates 
between  these  two  systems  and  results  in  an  unstable  evolution. 


To  begin,  suppose  at  first  that  c,  - c2  = 0 . 


For  any  initial  condition 


(.*,  (/„),;(:,(*„)),  the  total  mechanical  energy  is  conserved  in  Eq.  (5.2.15). 


A,  (t){mxX\  ( t ) + A:,*,  ( t ))  = 0 

d ( 1 1 ^ 

— -mxx2l{t)  + -kyxx(t) 
at  v 2 2 J 


= 0 


(5.2.18) 


jmlx?(t)  + ^klxf(t)  = Eto 


In  other  words,  the  trajectory  of  Eq.  (5.2.15)  is  an  ellipse  in  phase  space  parameterized  by 
the  initial  condition, 


yf(0  , A (0  x 


*.  1 

( \ 
mx 

KJ 

l2EJ 

(5.2.19) 


where  y,  (t)  = (t),  y2  (t)  = A,  (t).  A family  of  these  ellipses  is  shown  in  Fig.  5.1 

Similarly,  the  trajectory  of  Eq.  (5.2.16)  forms  an  ellipse  in  phase  space  given  by 


y,2(0  , y22  (0 


( \ 
m2 

l2^o  J 

2 E, 

v W 

= i 


(5.2.20) 


where  y,(t)  = x2(t),  y2(t)  = x2(t) . This  family  of  ellipses  is  shown  in  Fig.  5.2.  For 


purposes  of  illustration,  the  system  parameters  have  been  defined  as 
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k,>m,  for  Eq.  (5.2.15) 
k2<m2  for  Eq.  (5.2.16) 


( 


Figure  5.2  Stable  Phase  Plane  Trajectory  (2) 
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The  switching  rule  is  defined  via  Fig.  5.3.  Figure  5.4  depicts  a typical  trajectory  of  the 
switched  system. 


Figure  5.3  Switching  Rules:  Unstable  Case 


00 


Figure  5.4  Unstable  Phase  Plane  Trajectory 


Thus,  a pair  of  neutrally  stable  systems  can  generate  an  unbounded  trajectory  for  a 
particular  switched  system.  It  is  rather  simple  to  modify  this  argument  for  some  lightly 
damped  systems.  If  c,  * 0,  c2  * 0,  but  c,  « 1,  c2  « 1,  each  of  the  individual  systems 
is  asymptotically  stable  toward  the  origin.  With  the  same  switching  law,  the  phase  space 
diagram  for  the  system  in  Eq.  (5.2.17)  is  shown  in  Fig.  5.5. 
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00 


Figure  5.5  Unstable  Phase  Plane  Trajectory  with  Damping 


5.2.3  Stability  Conditions 


In  this  section,  it  is  shown  that  the  switched  system  used  to  approximate  the 
original  dynamical  system  does  not  exhibit  pathological  responses  of  the  type 
demonstrated  in  the  above  example.  Recall  from  Eqs.  (5.2.3)  and  (5.2.4)  that  only 
second  order  systems  of  ordinary  differential  equations  (ODE’s) 

[rri\x(t)  + [c]i(t)  + [&].*(/)  = 0 (5.2.21) 

that  can  be  written  as 


z(0  = [AM  t)  (5.2.22) 

in  the  event  that  the  forcing  function  is  identically  zero  are  considered.  It  is  not  difficult 
to  show  that  the  switched  system  described  by  Eqs.  (5.2.6)  and  (5.2.7)  is  equivalent  to 


k «>  - k lAh  (o 

h =s(ik_x,zikx(t)) 


(5.2.23) 


where  [pj  = 


This  leads  to  the  following  theorem, 
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Theorem:  For  any  switching  function  s ( • , • ) in  Eq.  (5.2.23)  and  initial  condition 
(i0,z0) , the  switched  system  in  Eq.  (5.2.23)  is  stable  provided  that 


M]r[^][0]+[0][^M^O  Vi 


(5.2.24) 


Proof:  The  proof  is  relatively  straightforward.  By  construction,  the  function 


V(z(t))  = V 


\ 


2l  (0 


= ^-{Zir(0  z\{t)} 


'[*]  0 ' 
0 [m] 


2i(0 

z2(t) 


(5.2.25) 


is  a Lyapunov  function  for  the  original  system  in  Eq.  (5.2.22). 


^m)  - \ 


f 


. T 


V 


- r'im]z,(n 


\ 

[h(t)[k]  l2(t)[m]]z(t) 

J 


(5.2.26) 


and  simplified  in  the  following  manner 


1 T 

—z(t) 

2" 


0 I 

i-l  r ;_i  r i-l 


[m]-'[k]  -[/«]"  [c] 


+ 


mm]} 


[*k,co 

[m]z2(t) 

0 / 
-[mT'[k]  -[m]-'[c] 


2i(0 

z2(0 


(5.2.27) 


=i{-zj(o[*]  (sro)[*i 

= ^(-z2r  (Ot^ki  (0  + (0[*k2  (0  _ ?2  (0[c]z2  (0) 

j?.(0 

jz2(0 

= ~2r(0[Ck(0 


= -^-{z,r(0  z2(0} 


o o 

o [c] 


(5.2.28) 

(5.2.29) 


(5.2.30) 


This  identity  can  be  rearranged  as  follows 
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d 

dt 


(k(zw»  - + 1 (.mu o) 

= i(zr(OMr[ek«)  + zr(0[GlMk(0) 

=^r<o{[/<]r[e]+[G][^]}z(o 
= -Lr«)[Ck<<) 


(5.2.31) 


where 


[0]  = 


'[*]  0 ' 

0 [m] 


and  [q  = 


0 O' 

0 [c] 


(5.2.32) 


Clearly,  V(zi(t))>0 


V{zt  (0)  = (0[*k  w (0  + (0[w]z2il  (0  > 0 


(5.2.33) 


and  moreover, 


dt 


(r{z,  (/»)  = Ui!  (OIGk,  (0  + 2.r  (Olfik,  (0) 


(5.2.34) 


= |(2,r(0[^]r[^][ek,(0+2,r(0[G][^][^k,(')) 

=^,r(o(Mfw][e]+[e]WM)2.v>  <5-2-35> 

= ^2,r(')[crk;(0 

However,  by  hypothesis  we  have 

M]r[^][0]  + [0][^M^  0 (5.2.36) 

This  implies  that  V is  a Lyapunov  function  for  all  i = 1...N  and  stability  follows  from 


Ref.  [12]  or  [54]. 
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5.3  Summary 

In  this  chapter  a modal  scheduling  technique  based  on  methodologies  for 
obtaining  reduced-order  models  via  switched  systems  was  derived.  While  switched 
system  modeling  has  provided  accurate  representations  in  the  control  of  various  nonlinear 
systems,  their  use  in  structural  dynamic  applications  has  been  limited.  This  chapter  has 
introduced  a reduced  model  definition  and  corresponding  logic  for  switching.  The 
models  used  in  the  switching  were  chosen  based  on  the  importance  of  the  modes.  The 
modal  participation  was  calculated  based  on  (1)  the  forcing  function  and  (2)  the  initial 
conditions  at  each  switch  time.  The  importance  of  stability  for  switched  systems  was  also 
discussed  and  a simple  example  of  an  unstable  switched  system  was  provided.  Stability 
criteria  for  the  proposed  switched  system  was  then  derived. 


CHAPTER  6 

SWITCHED  SYSTEMS:  NUMERICAL  EXAMPLES 


6.1  Introduction 

In  the  previous  chapter  a switched  system  model  and  corresponding  switching 
logic  was  derived  for  linear  structural  systems.  In  this  chapter  numerical  examples  are 
provided  to  study  and  demonstrate  the  usefulness  of  these  reduced  switched  system 
models.  The  first  numerical  example  considered  is  a simple  two-dimensional  decoupled 
spring/damper  system.  The  purpose  of  this  example  is  to  study  some  of  the  more 
fundamental  concepts  of  mode  switching  such  as  the  effects  of  reducing  switch  time  steps 
and  the  use  of  redundant  library  bases.  A second  and  more  complex  example  involves  a 
cantilevered  beam  system.  This  example  is  intended  to  demonstrate  the  potential  of 
mode  switching  in  linear  structural  analysis. 

6.2  Numerical  Example:  2-D  Spring  Damper  System 
6.2.1  Simulation  A 

The  system  considered  for  the  numerical  examples  in  this  section  is  a two- 
dimensional  spring  damper  system.  Since  the  system  is  defined  in  two  dimensions,  a 
reduced  model  only  consists  of  one  mode.  Thus,  the  modal  selection  method  described  in 
Section  5.2.1  is  simplified  considerably.  This  example,  however,  is  investigated  to  bring 
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to  light  aspects  of  switched  dynamical  systems  and  the  proposed  reduction  technique  that 
may  be  difficult  to  analyze  in  more  complicated  examples. 


Consider  the  dynamical  system  shown  in  Fig.  6.1. 
displacements,  the  decoupled  equations  of  motion  are  given  as 


m 

0 


0 

m 


k 0 jx 

0 *_W 


Figure  6.1  2D  Spring  Damper  System 


For  all  examples  in  Section  6.2,  the  system  parameters  are  defined  as  m 


lb/ft  and  c = 1 5 lb/ft/sec.  The  equations  in  state  space  form  are 
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Z{t)  = [A\  Z(t ) + F(t ) 


Assuming  small 


(6.2.1) 


► x 


= 5 slugs,  k = 500 


(6.2.2) 


(6.2.3) 


The  approximation  of  the  system  states  is  given  as 
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where  i again  refers  to  the  ith  interval  or  model.  In  this  first  simulation,  the  forcing 
function  is  defined  to  be  a sine  function  that  switches  between  the  positive  and  negative 
45°  angle  with  respect  to  the  X axis  as  shown  below. 


F(t ) = lOsin(lO^t) 


1/V21 

I/V2J 


F(t ) = lOsin(lO^r) 


t<  2 
2 < t < 6 
t>  6 


(6.2.5) 


For  this  demonstration,  four  basis  vectors  are  considered  in  the  basis  library  T. 


The  vectors  are 


2 

'll  42 

3 

' 1/V2  ’ 

4 

' 1/2  ' 

j 

’V3/2" 

¥ = 

_1/V2_ 

¥ = 

_-l  / V2_ 

¥ = 

_V3/2_ 

¥ = 

_ — 1 / 2 _ 

(6.2.6) 


and  correspond  to  the  45°,  -45°,  60°,  and  -30°  directionals  with  respect  to  the  X axis, 
respectively.  Note  that  the  presuperscripts  identify  the  modal  vector  in  F.  For  this 
example,  the  presuperscript  1 is  reserved  for  the  full-order  space  or  the  identity  matrix. 
The  vector  spaces  defined  in  Eq.  (6.2.6)  and  the  forcing  function  spaces  are  shown  in  Fig. 
6.2.  Note,  also,  that  the  forcing  function  lies  in  2y/  or  3 y/  depending  on  the  time.  One 
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would  expect  that  for  '/'defined  as 

T = [2y/  V V V]  (6.2.7) 

the  model  switching  logic  would  use  2 y/  for  the  first  six  seconds  of  the  simulation  and 

use  mode  3 y/  for  the  remaining  time  in  the  simulation. 

Simulation  results  for  the  full-order  model  or  the  “truth”  model  are  shown  in  Figs. 
6.3  and  6.4.  These  consist  of  X position  and  velocity  data.  Excited  motions  are  shown 
for  the  time  interval  between  zero  and  two  seconds,  they  dampen  out  in  the  interval 
between  two  and  six  seconds  and  restart  after  six  seconds. 


Figure  6.3  Full-Order  Model  - X Displacement 


Figure  6.4  Full-Order  Model  - X Velocity 
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In  the  reduced  simulations,  the  switch  interval  is  set  to  0.05  seconds  and  E.  is  set 

7 nun 

to  0.9.  The  modal  selection  algorithm  calculates  an  error  indicator  value,  E{  for  each 
modal  choice  and  the  mode  with  the  largest  value  is  chosen  as  the  mode  basis  in  the 
simulation  for  that  interval.  This  assumes  that  the  corresponding  Et  is  greater  than  Emin.  If 
this  is  not  true,  then  the  full-order  model  is  used. 

Full  and  reduced  simulations  are  compared  in  Figs.  6.5  and  6.6.  For  the  first  six 

seconds  the  second  mode  basis  2 y/  is  selected  and  the  reduced  system  matches  exactly 


Figure  6.5  Full  & Reduced-Order  Model  [2  3 4 5]  - X Displacement 


Figure  6.6  Full  & Reduced-Order  Model  [2  3 4 5]  - X Velocity 
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with  the  full  system  until  a time  of  six  seconds.  At  the  sixth  second,  an  error  initiates  and 
appears  to  dampen  out  after  about  the  eighth  second.  The  error  is  initiated  because  the 
forcing  function  (sine  function)  is  zero  at  six  seconds,  leaving  the  initial  conditions  as  the 
only  source  to  base  the  mode  selection  on.  At  the  sixth  second  the  minute  initial 
conditions  for  that  interval  reside  in  the  orthogonal  space.  Thus,  the  algorithm  selects  the 
second  basis  until  the  next  interval  at  6.05  seconds  is  evaluated.  At  this  time,  the  forcing 
function  is  no  longer  zero  and  resides  in  the  third  basis,  so  the  algorithm  then  chooses  the 
third  basis  based  on  the  evaluation  of  the  corresponding  Er 

Figures  6.7  and  6.8  show  semilog  plots  of  the  error  indicator  value  for  each  of  the 
four  modal  choices.  Since  modes  two  and  three  form  an  orthonormal  basis,  it  is  easy  to 
see  that  their  values  for  Ei  should  sum  to  one.  This  is  true  for  modes  four  and  five  as 
well.  In  Fig.  6.7,  Et  for  mode  two  shows  a residual  component  when  the  forcing  function 
is  switched.  This  is  due  to  the  minute  initial  conditions  that  reside  in  the  space  of  mode 


two  and  could  not  dampen  out  when  the  basis  is  switched  to  mode  three. 


Figure  6.7  Error  Indicator  for  Modes  2 and  3 
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Figure  6.8  Error  Indicator  for  Modes  4 and  5 


6.2.2  Simulation  B 


In  this  next  simulation,  the  third  mode  is  removed  from  T , resulting  in  the 
following  basis  library 

*¥  = [2  y/  4y/  5 y/~\  (6.2.8) 

With  this  basis  library,  it  is  again  expected  that  the  first  six  seconds  of  the  simulation 
would  use  2 \j/  and  again  be  exact.  Depending  on  the  minimal  accepted  error  indicator 

value,  the  algorithm  would  select  either  5y/  or  default  to  the  full-order  model  for  the 
remainder  of  the  simulation.  For  demonstration  purposes,  Emin  is  defined  low  enough  so 

that  s\j/  is  utilized.  Fig.  6.9  shows  the  results  of  the  X displacement  for  the  full  and 

reduced  models.  As  expected,  the  reduced  simulation  uses  3 y as  the  basis  after  six 
seconds  and  the  corresponding  results  are  not  very  accurate. 

In  an  attempt  to  study  a particular  consequence  of  the  model  switching  concepts, 
another  simulation  is  conducted  to  compare  output  that  results  for  smaller  switch 
intervals.  Recall  in  the  previous  simulation  that  a switch  interval  of  0.05  was  used.  In 


77 


the  next  simulation,  the  same  system  parameters  are  used  except  that  the  switch  interval 
is  reduced  from  0.05  to  0.01  seconds.  Figure  6.10  shows  the  decrease  in  error  when  the 
switch  time  is  reduced.  The  decrease  in  error  is  attributed  to  the  reduced  amount  of 

simulation  time  spent  in  2 y/  at  the  interval  starting  at  six  seconds. 


x 10'3 


Figure  6.9  Full  & Reduced-Order  Model  [2  4 5]  - X Displacement 


x 10'3 


Figure  6.10  X-Displacement  Error  - Switch  Interval  = 0.05  & 0.01 
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6.2.3  Simulation  C 


Another  consequence  of  the  switched  system  formulation  is  the  error  that  can  be 
induced  when  the  initial  conditions  of  the  interval  are  large  and  the  basis  that  is  chosen 
for  that  interval  is  orthogonal  or  almost  orthogonal  to  the  initial  conditions.  (Recall  Eq. 
(5.2.11))  If  a basis  is  chosen  that  does  not  span  the  space  of  the  initial  condition,  then 
there  is  no  means  for  that  initial  condition  to  dynamically  respond.  This  can  be  seen  in 
Eq.  (5.2.6).  In  the  previous  discussions,  the  transients  were  almost  completely  dissipated 
before  the  forcing  function  was  switched.  This  means  that  the  residual  error  at  the  time  of 
the  switch  was  quite  low.  In  the  next  set  of  simulations,  the  time  the  forcing  function  is 
switched  is  set  to  be  before  the  transients  have  had  time  to  dissipate  significantly.  At  the 
interval  when  the  forcing  function  is  switched,  Rt(t)  is  still  dominated  by  the  forcing 

function.  The  model  is  therefore  chosen  based  on  the  forcing  function,  but  the  initial 
condition  for  the  interval  is  significant.  Since  the  system  dynamics  occurs  in  the  space  of 
the  forcing  function,  any  portion  of  the  initial  condition  not  spanned  by  the  model  can  be 
thought  of  as  a residual  error. 

Figure  6.11  shows  the  full  and  reduced-order  models  simulations  with  forcing 
function  defined  as 


F(t)  = lOsin(lCto) 


fi/Vil 

I1/V2J 


t <2 


F(t)  = lOsin(lO^t) 


(6.2.9) 


t >2.5 
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Note  that  the  forcing  function  is  only  zero  for  2<t<2.5  as  opposed  to  the  previous 
simulation  where  it  was  zero  for  2<t<6.  The  switching  interval  is  again  0.05  seconds  and 
Emin  is  set  at  0.75.  The  basis  library  set  is  given  as 

= [*  y/  4y/  5y/]  (6.2.10) 

For  this  example,  the  switching  logic  selects  5 y/  for  the  first  two  seconds  and 

expectedly,  the  reduced  simulation  is  slightly  in  error.  For  t> 2.5  seconds,  y/  is  selected 
as  the  basis  and  the  residual  error  is  easily  seen.  This  is  because  the  space  of  the  model 
3 y/  no  longer  completely  spans  the  space  of  the  initial  conditions.  It  should  also  be  noted 
that  the  residual  error  shown  in  Fig.  6.1 1 is  extremely  large  because  of  the  nature  of  the 
two-dimensional  problem.  This  example  is  given  to  simply  demonstrate  a consequence 
of  the  approximation  defined  in  Eq.  (5.2.6).  As  systems  get  more  complex  and  the  model 
bases  are  depicted  by  more  than  one  mode,  these  errors  become  less  gross. 


x 10’3 


Figure  6.11  Switch  Intervals  with  Large  I.C.s  - X Displacement 
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Figure  6.12  shows  the  effect  of  reducing  the  switch  interval  and  increasing  Emin  on 
the  residual  error.  These  simulations  were  conducted  using  the  same  forcing  function  and 
library  basis  as  the  previous  example.  Intuitively,  one  might  expect  that  as  the  switch 
interval  is  decreased  and  as  the  error  tolerance  is  increased,  the  residual  error  should 
decrease.  However,  this  is  not  always  true.  As  the  switch  interval  is  decreased,  the  initial 
condition  of  interval  associated  with  the  start  of  the  second  forcing  function  may  be 
altered,  but  does  not  necessarily  move  toward  zero.  Similarly,  increasing  Emin  will  only 

improve  the  dynamic  space  and  does  not  affect  or  reduce  the  space  contained  by  the 
initial  condition. 


x 10‘3 


Figure  6.12  X Displacement  Errors  for  Switch  Interval  with  Large  I.C.s 


6.2.4  Simulation  D-l 

The  next  variation  of  this  first  numerical  example  addresses  issues  associated  with 
model  or  basis  redundancy.  That  is,  what  happens  when  the  vectors  in  the  basis  library 
overlap  or  span  similar  spaces.  Previously,  the  basis  library  was  defined  by  arbitrarily 
selecting  unit  vectors  with  particular  angles.  Now,  the  contents  and  size  of  the  basis 
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library  are  defined  by  a new  parameter,  termed  the  space  density  (SD).  This  literally 
means,  how  dense  or  how  heavily  populated  'F  is.  In  this  example,  SD  0 refers  to 
*F  containing  only  the  X and  Y axis.  SD  1 refers  to  a basis  that  contains  SD  0 plus  the 
spaces  or  angles  that  bisect  the  X and  Y axis.  As  shown  in  Fig.  6.13,  these  spaces  are 
simply  the  45°  unit  vectors  above  and  below  the  X axis.  Figure  6.14  shows  SD  2 which 
contains  SD  1 plus  the  unit  vectors  that  bisect  the  spaces  in  SD  1 . This  means  that  the 
angle  between  existing  spaces  is  22.5°  for  SD  2.  As  the  SD  increases  from  0, 1,  2 etc.,  the 


1 

0.8 
0.6 
0.4 

| 0.2 

1 

b 0 
I 

>- 

-0.2 
-0.4 
-0.6 
-0.8 

0 0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1 

X-  Direction 

Figure  6.13  Space  Density  - 1 
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Figure  6.14  Space  Density  - 2 
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angle  between  the  basis  vectors  decreases  from  90°,  45°,  22.5°  etc.,  and  the  number  of 
spaces  in  1 ¥ increases  from  2,  4,  8 etc.  Figure  6.15  shows  the  heavily  populated  space  of 
SD  7,  where  the  angle  between  the  spaces  is  0.7031°. 
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Figure  6.15  Space  Density  - 7 

In  the  next  simulations,  the  forcing  function  switches  are  again  separated  by  0.5 


-0.8 


0.4  0.5  0.6 

X-  Direction 


seconds,  so  that  the  residual  error  induced  by  the  initial  conditions  is  again  accentuated. 
The  forcing  function  is  defined  as 


F(t ) = lOsin(lO^t) 


f 0.93971 
{o.3420j 


F(t ) = lOsin(lO^t) 


| 0.9397  1 
{-0.3420  { 


t < 2 

2 < t < 2.5 
r >2.5 


(6.2.11) 


where  the  vector  directions  correspond  to  +/-  20°  about  the  X axis.  The  switch  interval  is 
given  as  0.005  seconds  and  Emin  is  set  at  0.90.  Figure  6.16  shows  the  velocity  error  in  the 
X direction  for  SD’s  2-5  and  Fig.  6.17  shows  the  same  error  for  SD’s  5-8.  A large  jump 
in  error  is  seen  at  2.5  seconds  when  the  forcing  space  or  direction  is  switched.  The  jump 
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Figure  6.16  Velocity  Errors  for  Space  Densities  2-5 
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Figure  6.17  Velocity  Errors  for  Space  Densities  5-8 


is  essentially  due  to  the  incorrect  basis  used  at  the  interval  associated  with  2.500<?<2.505. 


For  this  interval,  the  simulation  logic  chooses  the  basis  associated  with  the  +20°  forcing 


function  as  opposed  to  the  -20°  forcing  function.  Logic  can  be  added  to  reduce  this  error, 


however,  the  primary  scope  of  this  work  is  to  identify  consequences  of  switched  systems. 


Corrections  for  negative  consequences  may  be  addressed  at  a later  time. 


Intuitively,  as  the  library  basis  space  becomes  more  dense,  the  results  should 


become  more  accurate.  Figures  6.16  and  6.17  show  that  this  is  partially  or  globally  true. 
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In  Fig.  6.16  SD  3 does  not  appear  to  be  any  better  than  SD  2 and  in  Fig.  6.17  SD  7 yields 
better  results  than  SD  8.  This  issue  will  be  discussed  in  Section  6.2.6. 

6.2.5  Simulation  D-2 

In  the  next  study,  the  effects  of  reduced  switching  times  with  models  that  have 
extremely  redundant  bases  is  investigated.  In  simulations  B and  C,  the  same  effect  was 
studied,  but  with  models  containing  only  a few  basis  vectors.  The  results  from  C implied 
that  if  I/7  is  a fairly  poor  set  of  basis  vectors,  and  Emin  is  set  low  enough,  reducing  the 

switching  time  only  allows  switching  between  a bad  set  of  basis  vectors  more  often  and 
does  not  necessarily  improve  the  simulation  errors.  Will  this  be  the  case  if  the  library 
basis  is  more  dense?  This  is  the  question  that  is  addressed  next. 

For  Figs.  6.18-6.22  the  system  parameters  are  identical  to  those  defined  in 
simulation  D-l  except  that  the  switch  interval  value  is  run  at  0.01,  0.005,  and  0.001 
seconds.  The  higher  SD’s  (5  and  6)  seem  to  follow  the  logic  that  as  the  switch  interval  is 
reduced,  the  error  decreases.  As  with  the  results  from  simulation  C,  the  lower  SD’s  (2-4) 
don’t  necessarily  follow  this  logic. 
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Figure  6.18  Space  Density  2 - Various  Switch  Steps 
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Figure  6.19  Space  Density  3 - Various  Switch  Steps 
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Figure  6.20  Space  Density  4 - Various  Switch  Steps 
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Figure  6.21  Space  Density  5 - Various  Switch  Steps 
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Figure  6.22  Space  Density  6 - Various  Switch  Steps 
6.2.6  Simulation  D-3 


In  the  previous  two  simulations  D-l  and  D-2,  consequences  of  increasing  the 
library  basis  density  and  reducing  the  switch  interval  were  discussed.  Both  of  these 
studies  contain  added  errors  (i.e.,  the  large  spikes  at  three  seconds)  due  to  the  inability  of 
the  logic  to  anticipate  a correct  model  switch  for  the  interval  associated  with 
2.500<?<2.505.  In  this  section,  this  error  is  removed  by  simply  zeroing  out  the  initial 
forcing  function.  With  this  done,  the  effect  of  increasing  the  basis  spatial  density  should 
be  more  clear.  Thus,  for  this  demonstration  the  forcing  function  is  given  as 


F(t ) = 10sin(10;h) 


0.9397  1 
-0.3 420 J 


t < 0.2 
t >0.2 


(6.2.12) 


and  the  simulation  is  only  run  for  3 seconds  instead  of  10.  The  switch  interval  is  defined 
to  be  0.01  and  E-  is  set  at  0.90. 
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Figures  6.23-6.25  show  the  velocity  error  in  the  X direction  as  SD  ranges  from  2 
through  10.  What  may  not  be  immediately  obvious  is  that  the  errors  in  Fig.  6.23  are  an 
order  of  magnitude  higher  that  those  in  Fig.  6.24.  Similarly,  the  errors  in  Fig.  6.24  are  an 
order  of  magnitude  higher  that  those  in  Fig.  6.25.  Consequently,  the  error  decreases  by 
an  order  of  magnitude  for  every  three  increases  of  the  spatial  density.  This  can  easily  be 
seen  in  the  semilog  plot  shown  in  Fig.  6.26.  Careful  examination  of  one  particular  group 
of  SD’s  (i.e.,  2,  3,  and  4 or  5,  6,  and  7)  shows  that  there  is  no  correlation  between 
increasing  the  SD  and  reducing  the  simulation  error  within  the  space  of  the  three  SD’s. 
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Figure  6.23  Velocity  Errors  for  Space  Densities  2-4 


x 10" 


Figure  6.24  Velocity  Errors  for  Space  Densities  5-7 
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Figure  6.25  Velocity  Errors  for  Space  Densities  8-10 
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Figure  6.26  Velocity  Errors  for  All  Space  Densities 


Explanation  of  these  somewhat  strange  results  requires  investigating  the  basis  selections 
for  all  the  SD’s. 

In  Fig.  6.27  the  basis  vectors  that  are  used  in  the  simulation  for  SD’s  2,  3,  and  4 
are  shown.  For  SD  2,  there  are  eight  models  or  bases  to  choose  from  and  the  simulation 
locks  onto  model  six.  Similarly,  for  SD  3,  there  are  16  models  to  choose  from  and  the 
simulation  locks  onto  model  1 1 and  for  SD  4,  there  are  32  models  to  choose  from  and  the 
simulation  starts  with  model  21  and  also  uses  model  22.  The  starting  models  for  SD’s  2, 
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3 and  4 (6,  11,  and  21,  respectively)  are  all  associated  with  the  space  defined  by  -22.5°. 
When  SD’s  5,  6 and  7 are  analyzed,  the  initial  models  chosen  (see  Fig.  6.28)  are 
associated  with  the  space  defined  by  -19.6875°.  The  corresponding  model  numbers  for 
SD’s  5,  6 and  7 are  40,  79,  and  157,  respectively.  Similarly,  when  SD’s  8,  9 and  10  are 
analyzed,  the  initial  models  chosen  are  associated  with  the  space  defined  by  -20.0391°. 
For  this  example,  (forcing  function  at  -20°)  every  fourth  SD  provides  a more  accurate 
basis  for  the  simulation  to  start  with. 


Figure  6.27  Basis  Selection  for  Space  Densities  2-4 


Figure  6.28  Basis  Selection  for  Space  Densities  5-8 
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6.3  Numerical  Example:  Cantilevered  Beam 


In  the  previous  section  simple  numerical  simulations  were  conducted  to  identify 
and  study  numerical  problems  associated  with  switched  system  modeling.  In  this  section, 
a numerical  example  is  used  to  demonstrate  the  potential  and  advantages  of  modeling 
linear  structural  systems  as  switched  systems.  Figure  6.29  shows  a cantilevered  beam 
which  is  modeled  using  10  elements  and  has  30  DOFs.  A forcing  function  is  applied  at 
the  outboard  node  and  varies  between  axial  and  transverse  directions.  In  the  examples 
provided  two  values  for  Emin  are  used,  0.9  and  0.9999.  The  switching  interval  is  chosen 


to  be  0.01  seconds. 
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Figure  6.29  Natural  Frequencies  of  the  Full-Order  System 


The  forcing  function  at  node  1 1 is  defined  in  a cyclic  manner  as  follows 


F(t)  = 10sin(10^t) 


= 10sin(10;tf) 


0 < t < 0.5 
0.5  < t < 1 

1 < t < 1.5 
1.5  < t < 2 


(6.3.1) 


and  repeats  after  the  initial  2 seconds.  For  simulation  purposes,  the  system  parameters 
are  defined  to  yield  low  system  frequencies  and  are  given  as  follows 
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1 = 1 

p = \ 000 


a = 0.01 
E = 1.5x1 06 


/ = 1x10' 


(6.3.2) 


where  / is  the  element  length  (ft),  a is  the  cross  sectional  area  (ft2),  I is  the  moment  of 
inertial  (ft4),  p is  the  density  (sl/ft3)  and  E is  the  elastic  modulus  (lb/ft2). 

Figures  6.30-6.35  summarize  simulation  results  associated  with  Emin  = 0.9. 
Figure  6.30  shows  the  magnitude  of  the  two  functions  [A]Z0i  and  F\t)  which  are  used  to 
determine  modal  importance.  During  the  period  when  the  force  is  applied  transversely, 
the  relative  importance  of  the  forcing  function  is  shown  to  be  greater  than  that  of  the 
initial  conditions. 
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Figure  6.30  Force  and  IC  Magnitudes 

Figure  6.31  shows  the  number  of  modes  used  for  the  simulation  as  a function  of 
time.  Recall  that  there  are  a maximum  of  30  modes  which  can  be  used.  Figure  6.32 
shows  the  importance  of  each  mode  in  gray  scale  as  a function  of  time.  In  this  figure,  the 
color  black  is  most  important,  and  white  is  least  important  importance.  The  modes  are 
ordered  such  that  the  transverse  modes  comprise  modes  1-20,  and  the  axial  modes 


comprise  modes  21-30.  As  expected,  the  transverse  modes  are  the  only  modes  selected  in 
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the  first  second.  During  this  period  motion  is  contained  solely  in  the  space  spanned  by 
the  transverse  modes.  At  the  beginning  of  the  2nd  second,  the  force  is  applied  in  the  axial 
direction.  It  is  clear  that  the  selected  modes  from  the  period  t=  1 to  t= 2 seconds  are  axial. 

Figures  6.33-6.35  show  the  X and  Y displacements  and  rotation  of  the  11th  node 
for  the  full  and  reduced  (switched)  system.  The  plots  demonstrate  that  with  Emin  set  to 
0.9,  some,  but  not  all  of  the  full-order  dynamic  characteristics  are  captured.  Recall  from 
Fig.  6.31  that  only  between  10  and  15  of  the  possible  30  modes  were  utilized,  so  some 
error  is  expected. 


Figure  6.31  Number  of  Modes  Used 


Mode 


Figure  6.32  Modal  Importance 
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Figure  6.33  X-Displacement  of  Node  11 


Figure  6.34  Y-Displacement  of  Node  1 1 


Figure  6.35  Rotation  of  Node  1 1 
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The  next  set  of  results,  Figs.  6.36-6.41,  is  associated  with  Emin  set  at  0.9999. 
Figure  6.36  shows  the  magnitude  of  the  two  functions  [A]Z()i  and  F(t)  and  Fig.  6.37 

shows  the  number  of  modes  used  at  each  switch  time.  As  expected,  the  number  of  modes 
used  is  increased  considerably.  During  the  first  second  when  the  force  is  applied  in  the 
transverse  direction,  nearly  all  of  the  twenty  modes  are  included.  After  the  first  second 
there  is  a slight  increase  in  the  number  of  modes  until  1.5  seconds.  Then,  nearly  all  the 
modes  are  required  in  the  simulations. 


Figure  6.36  Force  and  IC  Magnitudes 


Figure  6.37  Number  of  Modes  Used 
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Figure  6.38,  presents  the  modal  importance  in  gray  scale.  As  expected  the  results 
look  very  similar  to  those  for  the  lower  value  of  Emin . Figures  6.39-6.41  again  show  the 
X and  Y displacements  and  rotation  of  the  1 1 th  node  for  the  full  and  reduced  (switched) 
system.  For  the  higher  value  of  Emin,  however,  there  is  a drastic  improvement  in 
simulation  results  for  the  reduced,  or  switched,  system. 
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Figure  6.38  Modal  Importance 
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Figure  6.39  X-Displacement  of  Node  1 1 


96 


Figure  6.40  Y-Displacement  of  Node  11 


Figure  6.41  Rotation  of  Node  1 1 


6.4  Summary 

In  this  chapter  numerical  examples  were  provided  with  two  goals  in  mind:  (1) 
investigating  fundamental  concepts  of  mode  switching  and  (2)  demonstrating  the 
potential  of  modal  scheduling  for  linear  structural  systems.  A two-dimensional  spring 
damper  system  was  used  as  the  model  for  studying  the  fundamental  concepts  and  a linear 
cantilevered  beam  was  used  to  demonstrate  the  value  of  this  reduction  technique. 
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For  the  two-dimensional  spring  damper  system,  a number  of  simulations  were 
conducted.  Simulation  A showed  that  errors  can  be  induced  when  the  simulation  stays 
with  a particular  mode  basis  for  just  one  extra  switch  interval.  Simulation  B 
demonstrated  even  larger  errors  when  a particular  mode  basis  was  removed  from  the  basis 
library  and  Emin  was  set  low  enough.  This  allowed  the  selection  algorithm  to  use  a basis 
that  was  nearly  parallel  with  the  forcing  function.  A second  simulation  showed  that  this 
error  could  be  reduced  by  decreasing  the  switch  step  size.  Simulation  C showed  that 
large  errors  can  be  introduced  if  the  initial  conditions  at  the  ith  interval  are  large  enough 
and  the  selected  basis  is  almost  orthogonal  to  the  initial  conditions.  Simulations  D-l,  2, 
and  3 studied  the  effects  of  using  grossly  redundant  basis  libraries.  These  simulations 
showed  there  was  a global  trend  for  the  simulation  error  to  decrease  as  the  density  of  the 
basis  library  is  increased.  However,  on  a more  local  level,  there  did  not  seem  to  be  a 
trend.  For  this  particular  example,  every  third  increase  in  space  density  showed  an  order 
of  magnitude  improvement  in  velocity  error.  This  was  simply  because  every  third 
increase  in  space  density  provided  a better  starting  basis  for  the  simulation.  Note  that  this 
result  is  only  particular  to  this  example  and  the  corresponding  forcing  conditions. 

The  last  numerical  example  concentrated  on  showing  the  benefits  of  using  modal 
scheduling  and  switched  systems.  A cantilevered  beam  was  used  as  the  model  with  a 
forcing  function  applied  at  the  outboard  body  node.  The  forcing  function  was  prescribed 
to  move  between  the  transverse  and  axial  directions  of  the  beam.  Two  different  values 
for  Emin  were  used,  0.9  and  0.9999.  Plots  showing  the  modal  participation  of  each  mode 
as  a function  of  time  were  used  to  demonstrate  that  only  certain  modes  needed  to  be 
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retained  in  the  model  to  obtain  relatively  accurate  simulation  results.  For  Emin  - 0.9, 
between  five  and  20  of  the  possible  30  modes  were  used,  but  the  simulations  showed  that 
only  part  of  the  dynamics  were  captured.  With  the  higher  value  for  Emin , between  15  and 
30  modes  were  selected  and  the  results  were  much  more  accurate.  Clearly,  a tradeoff 
exists  between  simulation  accuracy  and  the  amount  of  reduction.  This  tradeoff  is 
parameterized  by  Emin . 


CHAPTER  7 
CONCLUSIONS 


This  study  investigated  aspects  of  the  dynamic  simulation  of  mechanical  systems 
with  the  intention  to  improve  on  current  associated  technologies.  Two  aspects  were 
considered:  1)  the  modeling  of  flexible  components  and  2)  the  stabilization  of  the 
constrained  governing  equations. 

For  nonlinear  multibody  systems,  the  flexibility  of  components  is  typically 
modeled  using  an  assumed  modes  approach.  Conventional  methods  often  utilize  mode 
vectors  defined  at  the  component  level.  This  work  described  an  alternative  method  that 
involved  defining  component  modes  at  a system  level.  These  P&A  methods  generate 
reduced  components,  that  when  assembled,  form  reduced  systems  that  preserve  specified 
modal  frequencies.  This  is  achieved  by  projecting  selected  system  level  normal  modes 
onto  the  components.  This  approach  was  extended  by  using  system  level  Krylov  modes 
instead  of  normal  modes.  The  Krylov  method,  demonstrated  the  preservation  of  low 
frequency  moments  and  generated  frequency  responses  similar  to  the  P&A  methods.  The 
advantage  to  the  Krylov  method  is  that  the  full-order  eigenvalue  problem  can  be  avoided 
in  the  generation  of  system  level  modes.  This  model  reduction  procedure  was  also 
enhanced  by  using  the  balanced  realization  theory  to  determine  the  relative  importance  of 
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each  mode.  System  level  modes  were  ranked  and  only  the  most  important  retained  before 
they  were  projected  onto  the  components. 

A completely  novel  approach  for  modeling  flexible  components  involving  the  use 
of  switched  systems  was  also  presented.  In  this  approach  linear  structural  systems  that 
were  subject  to  time  and  spatial  varying  forces  were  treated  as  switched  systems. 
Switching  logic  was  derived  that  was  based  on  mode  participation  factors.  At  each 
switch  interval  participation  factors  for  each  mode  in  a basis  library  were  calculated. 
Only  the  most  important  modes  were  retained  in  the  model  for  that  switch  interval.  The 
modal  participation  factor  represented  each  modes  orthogonality  with  respect  to  the  initial 
condition  and  forcing  function  at  the  beginning  at  each  switch  step.  A simple  testbed 
example  was  used  to  study  and  investigate  some  of  the  more  fundamental  issues  of 
modeling  linear  structures  as  switched  systems.  A simple  spring/damper  system 
demonstrated  a number  of  effects  including,  the  consequences  of  using  inadequate  library 
bases,  and  using  large  switch  intervals.  The  study  also  included  the  effects  of  using 
extremely  redundant  library  bases.  A more  notable  result  was  shown  in  the  example 
involving  a cantilevered  beam.  These  results  showed  that  the  switched  model  can  often 
be  reduced  when  the  forcing  function  lies  within  either  the  axial  or  transverse  space  of  the 
beam.  Results  also  demonstrated  that  with  a relatively  low  minimum  acceptable  value 
for  the  error  indicator,  Emin , a small  number  of  modes  were  included  in  each  model 
resulting  in  simulations  that  captured  only  portions  of  the  dynamic  response.  Increasing 
the  value  of  E ■ lead  to  the  use  of  more  modes  and  consequently  better  dynamic 


simulations. 
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This  work  also  studied  the  stabilization  of  constraint  violation  drift  for  dynamical 
systems  with  redundant  constraints.  The  resulting  governing  equations  for  these  systems 
are  of  the  differential-algebraic  form  and  their  solutions  are  subject  to  constraint  violation 
drift.  In  this  dissertation,  it  is  shown  that  the  violation  in  the  DAEs  can  be  cast  into  a 
nonlinear  control  form  using  Lie  algebra.  Within  this  framework  the  stabilization  can  be 
interpreted  as  a controller  that  seeks  to  achieve  asymptotic  output  tracking.  Stabilization 
techniques  based  on  penalty  methods  were  also  derived  by  utilizing  energy  integrals 
associated  with  the  nonlinear  system.  These  methods  also  lead  to  a framework  that  was 
amenable  to  new  stabilization  techniques.  Numerical  examples  involving  a HMMWV 
and  weapon  system  were  used  to  demonstrate  the  penalty  methods  as  well  as  an 
alternative  method  based  on  the  sliding  mode  control  law  and  the  penalty  methods. 

There  are  several  subsequent  extensions  of  the  work  presented  in  this  dissertation. 
The  library  bases  used  in  the  mode  switching  contained  normal  modes  of  vibration.  An 
alternative  approach  would  be  to  use  Krylov  modes  or  combinations  of  the  Krylov  and 
normal  modes  in  the  library  bases.  Corresponding  switching  logic  would  also  need  to  be 
developed.  A second  and  more  difficult  extension  of  mode  switching  is  to  consider 
nonlinear  systems.  This  would  require  that  new  modal  participation  factors  be  developed 
based  on  external  forces  as  well  as  constraint  and  inertial  forces.  Further  extensions  can 
also  be  made  in  the  area  of  constraint  stabilization.  Both  of  the  aforementioned 
frameworks  allow  for  new  families  of  stabilization  methods  to  be  developed.  Numerical 
analysis  these  stabilization  families  on  complex  dynamical  systems  can  now  be  easily 


performed  using  DADS. 
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